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Abstract 


Title  of  Dissertation:  MODELING  AND  CONTROL  OF 

DYNAMICAL  EFFECTS  DUE  TO 
IMPACT  ON  FLEXIBLE  STRUCTURES 

Qifeng  Wei,  Doctor  of  Philosophy,  1994 

Dissertation  directed  by:  Associate  Professor  W.  P.  Dayawansa 

Department  of  Electrical  Engineering 


In  the  first  part  of  this  dissertation,  we  consider  modeling  and  approxima¬ 
tion  of  impact  dynamics  on  flexible  structures.  A  nonlinear  model  is  developed 
through  Hertz  law  of  impact  in  conjunction  with  the  dynamic  equation  of  the 
flexible  structure.  We  have  analyzed  this  nonlinear  model  and  established  the 
existence  and  uniqueness  of  solutions  of  the  nonlinear  equation.  A  numerical 
method  is  developed  based  on  the  contraction  mapping  principle.  By  utilizing 
the  fact  that  impact  interval  is  very  short  in  general,  one  may  approximate  the 
transfer  functions  of  the  systems  to  which  the  impacting  bodies  belong  by  Taylor 
polynomials  of  low  order.  We  have  developed  the  first  and  second  order  approx¬ 
imations.  The  first  order  approximation  yields  a  special  function  which  can  be 
used  for  analytical  and  computational  purposes.  The  second  order  approxima¬ 
tion  leads  to  a  two-parameter  family  of  ordinary  differential  equations  of  which 
the  solutions  exhibit  universal  features  of  impact  problems.  Simulation  results  of 
various  examples  have  demonstrated  the  usefulness  of  the  developed  numerical 
method  and  approximation  methods. 


The  second  part  of  this  dissertation  is  devoted  to  control  of  impact  dynamics. 
We  have  formulated  and  studied  a  control  problem  where  a  linear  system  is  sub¬ 
jected  to  a  series  of  impact  forces.  The  impact  forces  are  treated  as  disturbances 
to  the  system  and  modeled  as  finite  duration  events  using  the  theory  developed 
in  part  one.  A  reasonable  control  objective  is  to  design  a  feedback  controller  to 
minimize  the  energy  transferred  from  the  disturbances  to  the  controlled  outputs 
in  the  L2  norm  sense.  Under  the  assumption  that  the  disturbance  information 
is  known  a  priori,  a  (sub) optimal  control  strategy  is  derived  based  on  dynamic 
game  theory.  We  have  shown  that,  by  taking  advantage  of  the  fact  that  the 
duration  of  each  impact  force  is  very  short  in  general,  we  can  derive  a  series  of 
approximate  solutions  of  the  nonlinear  problem.  The  higher  order  terms  may  be 
negligible  for  the  disturbance  attenuation  problem  in  some  applications.  Hence, 
the  approximation  with  the  leading  term  renders  a  linear  one.  A  (sub) optimal 
Hoe  controller  is  derived  and  a  procedure  to  compute  such  a  controller  is  given. 
The  (sub)  optimal  solution  is  naturally  associated  with  the  existence  of  a  sta¬ 
bilizing  periodic  solution  of  coupled  Riccati  equations.  Hamiltonian  theory  is 
employed  to  analyze  the  coupled  Riccati  equations.  Finally,  we  investigate  the 
digital  implementation  of  this  control  algorithm  by  using  a  sampled-data  con¬ 
troller.  We  have  shown  that  under  a  certain  sampling  condition,  the  controller 
structure  could  become  simpler  than  the  continuous  version. 
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Chapter  1 


INTRODUCTION 


This  dissertation  is  composed  of  two  parts.  The  first  part  (chapters  2,  3  and  4) 
is  concerned  with  modeling  and  approximation  of  impact  dynamics  on  flexible 
structures.  Impact  phenomena  have  interested  scientists  and  engineers  for  a  long 
time.  Numerous  attempts  have  been  made  to  model  accurately  the  dynamical 
effects  of  impact  between  two  or  more  objects  [1,  2],  Our  motivation  to  study 
impact  problems  comes  from  force  and  constrained  motion  control  in  various 
robotic  manipulations  [1-8].  One  popular  example  in  these  problems  is  to  use  a 
flexible  robotic  arm  to  catch  a  baseball,  play  tennis,  snatch  a  stationary  object 
from  a  table,  hammer  a  nail  into  a  board,  or  collect  space  debris  while  attached 
to  a  moving  space  vehicle.  A  physical  contact  between  the  manipulator  and 
the  object  must  occur  before  the  desired  force  and  motion  can  be  applied  to 
the  object.  To  model  the  effects  of  these  actions  on  a  flexible  arm,  the  impact 
dynamics  must  be  examined.  Salmatjidis  [3]  and  Chapnik  etc.  [4]  independently 
studied  impact  control  of  flexible  manipulators.  Both  of  them  applied  open-loop 
control  strategies  to  cancel  the  impact  forces  due  to  the  bandwidth  limitation  of 
actuators  and  sensors.  The  open-loop  control  scheme  requires  a  precise  model  of 
impact  force  in  order  to  achieve  a  good  performance.  Salmatjidis  used  the  stere- 
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omechanical  impact  model  [5]  by  assuming  that  the  impact  duration  is  negligible. 
Chapnik  etc.  used  an  impact  model  devised  by  Lee  etc.  [6],  which  accounts  for 
the  finite  duration  of  impact.  Experimental  results  from  both  studies  showed 
that,  during  the  contact  period,  although  it  is  very  short,  impact  dynamics  can 
be  a  significant  factor.  If  the  impact  phenomena  are  not  modeled  properly,  or  the 
manipulator  is  not  controlled  to  diminish  the  effect  of  impact,  and  result  could 
be  the  failure  of  the  manipulator  [7,  8,  9,  10].  Moreover,  since  high  precision 
control  of  robotic  manipulators  has  become  increasingly  important  in  a  variety  of 
industrial  applications,  e.g.  laser  beam  technology,  semiconductor  wafer  manu¬ 
facturing  etc.,  this  requires  paying  extra  attention  to  the  usual  dynamical  effects 
as  well  as  taking  into  consideration  otherwise  ignored  features  such  as  dynamical 
effects  due  to  impact. 

Consideration  of  displacement  and  use  of  Hertz  law  of  impact  at  the  region 
of  contact  seems  to  be  the  most  successful  approach  [11].  When  the  contact 
involves  a  flexible  structure  Hertz  law  of  impact  leads  to  a  nonlinear  integral 
equation  called  the  Hertz  equation,  which  incorporates  the  effects  of  local  elastic 
deformation  at  the  region  of  contact  [6].  This  model  has  been  widely  applied 
to  various  impact  situations  [5,  12,  13],  and  experimental  results  obtained  in 
[12,  13]  well  support  the  validation  of  this  equation.  However,  there  are  several 
drawbacks  associated  with  this  nonlinear  model. 

1)  Numerical  methods:  this  nonlinear  equation  does  not  admit  a  closed  form  so¬ 
lution  in  general,  hence  numerical  methods  must  be  developed  to  solve  this  equa¬ 
tion.  There  are  some  numerical  methods  such  as  Timoshenko’s  small-increment 
method  and  an  energy  method  devised  by  Zener,  Feshbach  and  Lee  etc.  More 
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efficient  methods  need  to  be  developed  due  to  the  popular  use  of  this  model. 

2)  Mathematical  analysis:  although  the  above  numerical  methods  are  success¬ 
fully  used  to  solve  the  nonlinear  equation,  these  methods  were  proposed  without 
establishing  their  convergence  or  even  existence  of  a  unique  solution.  Our  view¬ 
point  is  that  theoretical  analysis  is  necessary  for  both  proving  the  validation  of 
this  equation  and  developing  efficient  numerical  methods. 

3)  Computational  complexity:  in  all  existing  numerical  methods,  each  case,  e.g. 
varying  initial  velocities,  has  to  be  numerically  solved  separately.  Also,  a  fairly 
large  computational  burden  has  to  be  incurred  for  each  numerical  solution.  One 
of  the  goals  of  this  dissertation  is  to  attempt  to  solve  these  problems. 

Our  methodology  for  modeling  the  impact  dynamics  is  based  on  the  Hertz 
law  of  impact.  In  chapter  3,  we  derive  a  nonlinear  impact  model  through  Hertz 
law  of  impact  in  conjunction  with  the  dynamical  equation  of  a  flexible  structure 
which  involves  impact.  We  establish  the  existence  and  uniqueness  of  solutions  of 
the  Hertz  equation  by  applying  the  contraction  mapping  principle.  The  detailed 
proofs  of  both  local  and  global  solutions  will  be  given.  Chapter  4  addresses 
the  problem  of  numerically  solving  the  Hertz  equation.  Timoshenko’s  small- 
increment  method  [11]  yields  a  solution  within  any  degree  of  accuracy  provided 
that  sufficiently  small  time  steps  of  integration  are  used.  Hence  it  has  become 
the  basis  for  evaluating  other  approximation  methods.  Another  useful  method 
is  the  application  of  the  energy  method  devised  by  Zener  and  Feshbach  [14],  and 
applied  by  Lee  [6]  to  the  central  impact  of  a  sphere  on  a  simply  supported  beam. 
By  combining  the  principle  of  conservation  of  energy  during  the  impact,  and 
momentum  considerations  for  the  impacting  sphere,  the  resulting  velocity  of  the 
impacting  sphere  and  vibrational  energy  imparted  to  the  beam  are  determined. 
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Impact  force  was  assumed  to  be  a  sinusoidal  function  of  time.  This  method  is 
simple  and  fast,  but  not  reliable,  and  we  will  give  some  numerical  examples  to 
illustrate  this.  In  its  place  we  will  develop  a  numerical  method  using  successive 

Picard  iterations;  /o,  Pfo ,  PPfo, . ,  where  the  initial  condition  /0  is  obtained 

from  the  energy  method,  and  P  is  a  contraction  operator.  We  show  that  this 
method  is  faster  compared  to  the  small-increment  method,  more  accurate  than 
the  energy  method,  and  its  convergence  is  very  fast. 

Observe  from  numerical  solutions  that  the  impact  period  is  very  short  in 
general.  Therefore,  one  may  approximate  the  transfer  functions  of  the  systems 
to  which  the  impacting  bodies  belong  by  Taylor  polynomials  of  low  order.  We 
explicitly  carry  out  this  computation  in  the  cases  of  first  and  second  order  Taylor 
polynomial  approximations  of  the  transfer  functions.  We  show  that  in  the  case 
of  the  first  order  approximation  there  is  a  universal  ordinary  differential  equation 
that  describes  the  impact  behavior  completely.  Therefore,  one  can  numerically 
solve  this  equation  beforehand,  save  the  results,  and  use  it  to  predict  the  impact 
behavior  with  only  a  minimal  computational  burden.  In  the  case  of  the  sec¬ 
ond  order  approximation  there  is  a  two-parameter  family  of  ordinary  differential 
equations  that  govern  the  impact  behavior.  These  approximation  methods  are 
validated  via  numerical  examples. 

The  second  part  (chapter  5  and  6)  is  devoted  to  control  aspects  of  impact 
dynamics.  We  formulate  and  study  a  control  problem  where  a  linear  system 
is  subjected  to  a  series  of  impact  forces.  The  impact  forces  are  treated  as  the 
disturbances  to  the  system  and  modeled  as  finite  duration  events  using  the  theory 
developed  in  part  one.  Among  the  motivating  factors  is  the  need  to  study  the 
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control  problems  related  to  mechanical  systems  subject  to  impact  forces,  e.g. 
active  control  of  the  suspension  system  of  a  vehicle  against  irregularities  of  the 
road  [15,  16],  accurate  pointing  of  guns,  stabilization  of  an  antenna  on  the  space 
station  subject  to  impact  from  space  debris,  or  active  damping  of  vibrations 
of  flexible  structures  caused  by  impact  forces  [10,  17].  A  reasonable  control 
objective  in  all  these  problems  is  to  design  a  stabilizing  controller  to  minimize 
the  energy  transferred  from  the  disturbances  to  the  controlled  outputs  in  the 
L2  norm  sense.  This  problem  in  turn  can  be  studied  in  the  framework  of  H0 0 
control  theory. 

An  important  paradigm  in  control  synthesis  is  the  H0 0  control  problem  in¬ 
troduced  by  Zames  [18].  In  this  formulation,  the  disturbances  belong  to  a 
ball  in  a  certain  function  space,  and  a  quadratic  cost  function  is  minimized 
for  the  worst  disturbance  in  this  set.  Various  theories  of  H0 0  control  prob¬ 
lems  for  linear  systems  have  been  developed  so  far  by  many  researchers  (see 
e.g.  [19,  20,  21,  22,  23,  24]).  Solutions  to  these  problems  can  be  obtained  via 
frequency  domain  methods,  or  recently  developed  state-space  methods. 

Recently,  several  researchers  have  attempted  to  extend  control  results 
to  nonlinear  systems.  Ball  and  Helton  [25,  26],  from  a  viewpoint  of  operator 
theory,  have  discussed  control  theory  of  nonlinear  systems  for  the  first  time. 
Basar  and  Isidori  [27,  28]  study  the  connection  between  the  H0 0  control  theory  of 
nonlinear  systems  and  differential  game  theory.  In  this  setting,  one  is  naturally 
led  to  a  nonlinear  partial  differential  equation  known  as  the  Isaacs  equation. 
A  straightforward  application  of  the  theory  (see  e.g.,  [29])  to  the  case  of  a  plant 
modeled  by  an  affine  nonlinear  system  shows  that  once  a  solution  of  the  ap¬ 
propriate  Isaacs  equation  is  found,  a  (full-information)  feedback  law  providing 
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disturbance  attenuation  (in  the  L2  gain  sense)  can  be  computed  right  away.  Van 
der  Schaft  [30,  31]  analyzed  the  relation  of  the  L2  gain  between  nonlinear  sys¬ 
tems  and  their  linearization,  and  gave  a  sufficient  condition  for  the  existence 
of  smooth  Hoc  state  feedback.  In  addition,  using  the  dissipative  system  theory 
[32,  33],  Van  der  Schaft  has  discussed  the  relation  between  the  L2  gain  and  the 
Hamilton- Jacobi  equation  in  [34].  However,  his  approach  requires  the  equilib¬ 
rium  point  of  a  suitable  Hamilton  system  to  be  hyperbolic. 

The  control  of  nonlinear  systems  is  still  largely  open.  One  problem  is  the 
output  feedback  control.  Van  der  Schaft  derived  an  output  feedback  controller  for 
a  certain  class  of  nonlinear  systems  [35],  which  however  has  the  severe  drawback 
of  being  infinite-dimensional.  Much  more  effort  has  been  put  into  obtaining 
finite-dimensional  controllers  in  general  (see,  e.g.  [36]),  but  without  decisive 
answers  thus  far.  Another  main  issue  for  the  applicability  of  nonlinear  Hoc 
control  is  its  computational  complexity,  one  of  the  essential  requirements  in  all 
of  the  above  approaches  is  either  to  solve  a  Hamilton-Jacobi  equation  or  a  set  of 
completed  Hamilton-Jacobi  inequalities  [34,  35]. 

The  disturbance  attenuation  of  impact  forces  is  a  nonlinear  control  problem 
due  to  the  nonlinearity  of  the  impact  model,  which  is  different  from  the  prob¬ 
lems  discussed  in  [36,  34,  35].  It  is  not  affine  in  disturbance  input,  and  the 
linearization  around  the  equilibrium  does  not  exist.  Hence  this  nonlinear  prob¬ 
lem  cannot  be  solved  by  the  above  proposed  methods,  but  some  analysis  can 
be  carried  over  to  this  problem.  Our  approach  is  based  on  the  dynamic  game 
theory  [27,  29].  In  this  setting,  the  control  problem  naturally  becomes  a  mini¬ 
max  optimization  problem  of  a  cost  function  L(u,v)  with  two  players,  where  u 
is  the  control  and  v  is  the  disturbance.  Roughly  speaking,  the  player  u  tries  to 


6 


minimize  L(u,  v)  in  U,  while  the  player  v  tries  to  maximize  L(u,  v)  in  V  simul¬ 
taneously,  where  U,  V  are  the  constraint  sets  of  u  and  v  respectively.  [27,  29] 
show  that  the  (sub)  optimal  disturbance  attenuation  problem  has  a  solution  for 
a  given  7  >  0  if  the  minimax  optimization  problem  admits  a  saddle  point.  We 
show  that,  due  to  nature  of  impact  dynamics,  the  saddle  point  may  not  exist  in 
problems  involving  impacts.  Motivated  by  the  dynamic  game  approach,  if  the 
information  of  the  disturbance  v  is  assumed  to  be  known  a  priori  (e.g.  sensors 
can  predict  the  impact  velocity  before  impact),  one  can  seek  an  optimal  strategy 
by  using  this  information  such  that  a  certain  attenuation  level  7  is  achieved.  A 
procedure  is  given  to  compute  the  optimal  strategy  u(v),  and  the  optimal  atten¬ 
uation  level  7*.  Finally,  by  taking  advantage  of  the  fact  that  the  duration  of  each 
impact  force  is  very  short  in  general,  we  derive  a  series  of  approximate  models 
for  the  original  nonlinear  system.  In  many  cases  it  turns  out  that  the  higher 
order  terms  are  negligible.  Hence  a  special  linearization  of  the  nonlinear  system 
can  be  obtained  by  using  the  leading  term  as  an  approximation.  Thus,  in  these 
cases  the  nonlinear  problem  can  be  approximated  by  a  linear  one.  In  chapter  6, 
motivated  by  the  linear  impulsive  model  obtained  in  chapter  5,  we  formulate  a 
linear  control  problem.  In  the  infinite  horizon  case,  a  state-space  solution 
for  the  state-feedback  controller  design  problem  is  derived.  The  (sub) optimal 
JToo  controller  is  naturally  associated  with  the  existence  of  stabilizing  periodic 
solution  of  a  coupled  Riccati  equations  (one  is  a  differential  Riccati  equation, 
one  is  a  difference  Riccati  equation).  Hamiltonian  theory  is  employed  to  analyze 
the  coupled  Riccati  equations.  Finally,  we  investigate  the  digital  implementa¬ 
tion  of  this  control  algorithm  by  using  a  sampled-data  controller.  We  show  that 
under  a  certain  sampling  condition,  the  control  problem  can  be  converted  into 
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a  standard  discrete-time  H ^  problem.  Then,  the  state-space  solution  is  derived 
through  dynamic  game  theory.  The  structure  of  the  sampled-data  controller 
could  be  simpler  than  the  continuous  version.  An  impact  example  is  given  and 
discussed. 

In  the  final  chapter,  we  summarize  the  main  results  obtained  in  this  disser¬ 
tation,  and  identify  some  interesting  problems  for  future  research. 


8 


Chapter  2 


A  REVIEW  OF  PREVIOUS  WORK 


In  this  chapter,  we  will  review  some  previous  work  on  modeling  and  control  of 
impact  dynamics. 

2.1  Impact  Modeling 

Impact  phenomena  have  interested  scientists  and  engineers  for  a  long  time.  Nu¬ 
merous  attempts  have  been  made  to  accurately  model  dynamical  effects  of  im¬ 
pact  between  two  or  more  objects.  The  first  formulation  of  rigid-body  impact 
theory  is  due  to  Galileo,  which  the  focus  is  on  the  impulse-momentum  law  for 
rigid  bodies  and  involves  a  minimum  of  mathematical  difficulties.  For  perfectly 
elastic  impact  of  two  bodies,  the  law  of  conservation  of  energy  provides  the  sec¬ 
ond  relation  required  to  uniquely  determine  the  final  velocities  of  the  objects. 
Newton  later  introduced  a  coefficient  of  restitution  e  to  account  for  the  degree 
of  plasticity  of  the  collision  and  energy  loss  during  impact.  This  modified  im¬ 
pact  law  was  widely  used  in  many  impact  problems  [1,  2,  15,  16].  However,  it  is 
incapable  of  describing  the  transient  stresses,  forces,  or  deformations  produced, 
and  is  limited  to  a  specification  of  the  initial  and  terminal  velocity  states  of  the 
objects  and  the  applied  linear  or  angular  impulse.  The  theory  fails  to  account 
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for  local  deformations  at  the  contact  point  and  further  assumes  that  a  negligible 
fraction  of  the  initial  kinetic  energy  of  the  system  is  transformed  into  vibrational 
energy  of  the  colliding  bodies.  The  last  hypothesis  has  been  found  to  be  reason¬ 
ably  valid  for  the  collision  of  two  spheres,  but  not  for  collisions  involving  flexible 
structures  (beams,  rods,  plates  etc). 

The  first  satisfactory  analysis  of  the  stresses  caused  by  the  impact  of  two 
elastic  bodies  is  due  to  Hertz,  who  viewed  the  contact  of  two  bodies  as  an 
equivalent  problem  in  electrostatics.  A  solution  was  obtained  in  the  form  of  a 
potential  which  described  the  stresses  and  deformations  near  the  contact  point 
as  a  function  of  the  geometrical  and  elastic  properties  of  the  bodies. 

When  impact  involves  a  flexible  structure,  consideration  of  displacement  and 
use  of  Hertz’s  law  of  impact  at  the  region  of  contact  seems  to  be  the  most  suc¬ 
cessful  approach  [11],  Hertz  law  of  impact  leads  to  a  nonlinear  integral  equation 
called  the  Hertz  equation,  which  incorporates  the  effects  of  local  elastic  deforma¬ 
tion  at  the  region  of  contact  [6].  This  model  has  been  widely  applied  to  various 
impact  situations  [5,  12,  13],  and  the  experimental  results  obtained  in  [12,  13] 
well  support  the  validity  of  this  equation.  Unfortunately,  this  nonlinear  equa¬ 
tion  does  not  admit  a  closed  form  solution,  hence  numerical  methods  must  be 
developed  to  solve  this  equation.  Timoshenko  [11]  developed  the  so-called  small- 
increment  method  to  give  numerical  solutions,  and  it  has  become  the  basis  for 
evaluating  other  approximation  methods.  Some  other  approximation  methods 
also  give  very  satisfactory  results  [6,  38].  One  of  them  is  the  application  of  the 
energy  method  devised  by  Zener  and  Feshbach  [14],  and  applied  by  Lee  [6]  to 
central  impact  of  a  sphere  on  a  simply  supported  beam. 
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2.1.1  Stereomechanical  impact 

The  classical  theory  of  impact,  called  stereomechanics,  is  based  primarily  on  the 
impulse-momentum  law  for  rigid  bodies  and  involves  a  minimum  of  mathemati¬ 
cal  difficulties  in  its  formulation.  For  perfectly  elastic  impact  of  two  bodies,  the 
law  of  conservation  of  mechanical  energy  provides  the  second  relation  required 
to  uniquely  determine  the  final  velocities  of  the  objects.  When  the  impact  pro¬ 
duces  a  permanent  deformation  this  relation  is  replaced  by  the  introduction  of  a 
coefficient  of  restitution  e  for  the  process.  This  coefficient  purports  to  describe 
the  degree  of  plasticity  of  the  impact,  and  is  usually  defined  as  the  ratio  of  final 
to  initial  relative  velocity  components  of  the  striking  objects  in  the  direction 
normal  to  the  contact  surfaces.  This  impact  theory  has  been  used  in  various 
applications,  e.g  central  impact  of  two  rigid  bodies.  Assume  the  two  rigid  bod¬ 
ies  of  masses  mi  and  m2  respectively  make  initial  contact  with  each  other  at  a 
point  on  the  line  connecting  their  centers  of  gravity  with  velocities  vw  and  u2o 
respectively.  By  the  use  of  the  moment  conservation  law,  the  velocities  of  the 
two  bodies  after  impact  are  given  by  [5], 

(l  .  ^2(^10  —  ^20)  /r.  ,  -.N 

vij  =  v10-(l  +  e) - ; -  (2.1.1) 

mi+m2 

,  M  ,  iml(v10-V2o) 

v2f  —  v2o  +  (1  +  e) - - - 

m  1  +  m2 

—  V20  +  K2(v  10  —  V20),  (2.1.2) 

where  e  is  the  coefficient  of  restitution  and  assumed  to  be  a  constant.  For  a 
given  e  and  initial  velocities,  the  final  velocities  of  contacting  objects  can  be 
determined  immediately  from  these  equations. 

However,  there  are  some  limitations  associated  with  this  method,  especially  for 
impact  involving  a  flexible  structure.  It  is  incapable  of  describing  the  transient 
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stresses,  forces,  or  deformations  produced,  and  the  method  further  assumes  that 
a  negligible  fraction  of  the  initial  kinetic  energy  of  the  system  is  transformed 
into  vibrational  energy  of  the  colliding  bodies. 

2.1.2  Hertz  law  of  impact 

The  first  satisfactory  analysis  of  the  stresses  at  the  contact  of  two  elastic  bodies 
is  due  to  Hertz,  who  viewed  the  contact  of  two  bodies  as  an  equivalent  problem 
in  electrostatics.  A  solution  was  obtained  in  the  form  of  a  potential  which 
described  the  stresses  and  deformations  near  the  contact  point  as  a  function  of 
the  geometrical  and  elastic  properties  of  the  bodies.  The  Hertz  law  of  impact  is, 

<*(<)  =  (2.1.3) 

where  a(t)  is  the  relative  approach,  i.e.  the  difference  between  the  displacements 
of  the  bodies,  and  K  is  the  Hertz  constant  [5,  39], 

K=^{q/(Q1+Q2VA  +  B)},  (2.1.4) 

where  q,  A  and  B  are  constants  which  depend  on  the  local  geometry  of  the  region 
of  contact,  for  a  sphere  of  radius  Ri  and  plane  surface  contact  case,  A  —  B  =  —■ 
and  q  —  7T1/3,  and, 

Qi  =  (i  -  iA  Qi  =  (i  -  i4)/e2Ki  (2.1.5) 

where,  Hi  and  Ei  are  the  Poisson  ratio  and  Young’s  modulus  for  the  two  bodies 
respectively. 

This  impact  law,  although  both  static  and  elastic  in  nature,  has  been  widely 
applied  to  various  impact  situations  and  the  experimental  results  obtained  well 
support  the  validity  of  this  equation. 


12 


2.1.3  Central  impact  of  a  sphere  on  a  simply  supported 
beam 

Timoshenko  was  first  to  use  the  Hertz  law  of  impact  to  model  the  impact  dy¬ 
namics  on  a  flexible  beam.  He  formulated  and  solved  numerically  the  problem 
of  central  impact  of  a  sphere  on  a  simply  supported  beam  [11,  5].  The  impact 
problem  can  be  formulated  as  follows;  a  beam  is  struck  transversely  by  a  spher¬ 
ical  mass  m  with  initial  moving  velocity  Vo.  We  further  assume  that  the  Hertz 
law  of  impact  is  valid  i.e. 

Q  =  K[f(t)fl\  (2.1.6) 

where  a  is  the  relative  approach  of  the  striking  body,  f(t)  is  the  contact  force. 


Figure  2.1:  A  sketch  of  the  displacement 

The  relative  approach  is  the  difference  between  the  displacement  of  the  beam 
and  the  contacting  body,  measured  from  the  instant  of  initial  contact.  Hence 

a  =  s(t)  —  w(x*,t),  (2.1:7) 


13 


where  w(x*,t)  is  deflection  of  the  beam  at  the  point  of  contact  x*,  s(t)  is  the 
displacement  of  the  ball  under  of  the  contact  force  f(t),  and  is  given  by, 


(t)  =  v0t  -  —  [  f(r)(t  -  r)dr. 
m  Jo 


(2.1.8) 


From  the  equations  (2.1.6)-(2.1.8)  we  obtain  the  nonlinear  integral  equation, 

K[f(t)]2/3  —  V(jt  —  f  f(j)(t  -  r)dr  -  w(x*,  t).  (2.1.9) 

m  Jo 

The  forced  vibrations  produced  in  the  beam  by  the  varying  force  fir)  at  the 
point  of  contact  are  expressed  in  terms  of  the  normal  modes  of  vibration.  In 
the  case  of  symmetrical  vibration  of  a  simply  supported  beam,  the  normalized 
characteristic  deflections  [11]  are  \j2/lsin(mrx/l)  for  n  =  1,3,5,  etc.,  where  x 
is  the  coordinate  of  position  along  the  beam  of  length  l.  The  corresponding 
frequencies  of  the  natural  modes  of  vibration  are  wn ,  related  by  the  equation 
wn  =  n2Wi- 


The  central  deflection  w(x*,t)  due  to  the  forced  vibrations  is  given  by 

2  /•*  sinwn{t-T ) 


w(x 


ra=l,3,5 


-dr 


wn 


(2.1.10) 


where  q  is  the  constant  mass  per  unit  length  of  the  beam,  ql/ 2,  denoted  by  M, 
is  called  the  reduced  mass  and  in  this  particular  case  is  equal  to  half  the  mass  of 
the  beam.  Substituting  equation  (2.1.10)  in  equation  (2.1.9)  gives  the  integral 
equation  for  the  impact  force  from  which  the  motion  can  be  determined,  that  is 

K[f(t)]2/*=v0t-^-  j  f(r)(t-r)dr-  ~  f  /(T)sww"|* — ^dr.  (2.1.11) 

771  JO  n=l,3,5  ^  ^ n 

The  nonlinear  term  K[f(t)]2/3,  precludes  any  closed  form  solution  for  (2.1.11). 
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2.2  Numerical  Methods 


There  are  some  numerical  methods  available  in  the  existing  literatures  [6,  38,  40]. 
Timoshenko’s  small-increment  method  [11]  is  usually  reliable,  but  could  turn  out 
to  be  time  consuming.  Another  useful  method  is  the  application  of  the  energy 
method  devised  by  Zener  and  Feshbach  [14],  and  applied  by  Lee  [6]  to  central 
impact  of  a  sphere  on  a  simply  supported  beam.  This  method  is  simple  and  fast, 
but  not  reliable  as  we  will  illustrate  via  some  numerical  examples. 


2.2.1  Small-increment  method 


In  this  section  we  will  illustrate  Timoshenko’s  small-increment  method  (  see  [11] 
for  details).  Recall  that  the  Hertz  equation  is  given  by 


I  pi  OO  I  « 

K[f{t)]2/3  =  v0t--  I  f(r){t-T)dT-  E  T7L 

to  Jo  i=i>3)5  M  Jo 


)rfr(221) 
Wi 


Small-increment  method  divides  the  time  interval  [0,  t)  into  n  sub-intervals 
with  increments  At,  and  the  impact  force  is  assumed  to  be  constant  during 
each  sub-interval.  Thus  for  the  nth  time  interval  t  —  nr,  equation  (2.2.1)  can  be 
written  as 


K[fn]2/3  =  v0nAr  -  E  Dn-j+lfj  ~ 

m  j= i 

1  ^  t  ^  cosWi(n- j)Ar  -  cosWi(n- j +  l)Ar,n  n 

1V1  j= 1  i=l,3,5  Wi 

where  Ar  is  conveniently  chosen  as  some  small  fraction  of  the  fundamental 
period  of  vibration  of  the  beam.  The  quantity  Y2j=i  Dn-j+ifj  arises  from  the 
double  integration  term  and,  for  a  linear  continuous  approximation  of  the  force- 
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time  curve  may  be  expressed  as 

^  Ai-j+i/j  =  2[(n  —  1  )fi  +  (n  —  2 )(/2  —  /i)  4 - 

j=i 

(n  ~  i)(/j  ~  fj-i  +  /i-2  -  /j— 3  H - T  /i)  H - b 

(/n— 1  —  fn—2  +  fn- 3  —  •  •  •  ±  /i)]  + 

l/3(/n  -  fn—1  +  /„-2  -  '  *  •  ±  /l)  (2.2.3) 

The  accuracy  of  the  results  and  the  labor  of  computation  depend  upon  the 
proper  evaluation  of  the  term  coswAn~i)^T~^°swM~i+l)h-r  being  directly  pro¬ 
portional  to  the  number  of  modes  included  and  the  fineness  of  the  time  inter¬ 
val  chosen.  Procedural  techniques  for  the  step-wise  approximation  are  given 
in  reference  [41].  Although  this  method  is  very  lengthy,  it  is  considered  in  the 
literature  as  the  standard  numerical  method  for  solving  the  Hertz  equation. 

2.2.2  A  solution  by  approximating  impact  force 

Lorenertz  [42]  developed  a  concise  method  of  solution  by  assuming  that  the 
duration  of  impact  is  small  compared  with  the  period  of  the  fundamental  mode 
of  vibration  of  the  beam  and  also  that  only  fundamental  model  is  excited.  The 
first  assumption  is  justified  in  most  impact  problems  and,  moreover,  is  easily 
checked,  but  the  decision  where  or  not  it  is  permissible  to  neglect  the  effect  of  the 
higher  modes  of  vibration  of  the  beam  is  more  difficulty.  The  latter  assumption 
is  equivalent  to  replacing  the  beam  by  the  one  degree-freedom  system  of  a  mass 
on  a  spring  [6]. 

The  method  of  solution  is  to  assume  a  sinusoidal  contact  variation  /(f)  = 
Psin(Xt)  which  enables  the  integrals  in  equation  (2.2.1)  to  be  evaluated,  thus 
transforming  it  into  an  algebraic  equation.  Values  of  P  and  A  are  required  to 
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satisfy  this  equation  for  all  t  during  the  impact.  Since  f(t)  =  Psin(Xt)  cannot 
be  negative,  contact  occurs  for  0<  A t  <  7r,  so  that  the  duration  of  contact 
is  Td  =  7t/A (sec).  Substituting  the  sinusoidal  variation  for  f(t)  and  and 
evaluating  the  integral,  the  equation  (2.2.1)  becomes 


■  s  vm  P  sinXt .  P  Xsinw it  —  WisinXt 

K[Psmxtr  =  -  -(-x  -  — )  -  M  roi(V-W)  •  (2'2'4) 

Equation  (2.2.4)  is  still  too  complicated  for  direct  solution.  Some  simplifications 
can  be  made  by  writing  w\t  for  sinw\t ,  since  w\t  remains  small.  In  addition 
[sin2/3  At]  can  be  approximated  by  AsinXt\  where  A  is  chosen  so  that  the  new 
contact  force  variation,  which  makes  a  —  KP2/3AsinXt,  produces  the  same 
total  change  of  momentum  in  the  striking  mass  as  did  PsinXt.  This  condition 
expressed  mathematically  becomes 


/  sinOdO  —  A  f  sin2/30d0 
Jo  Jo 

which  gives  A  —  1.093  ,  inserting  these  simplifications,  Equation  (2.2.4)  be¬ 
comes 


KP2P  A  sinXt 


P  ,t  sinXt .  P  Xt  —  sinXt  ,  . 

Vot - (t  -  -7^)  ~  2.2.5 

m  X  X1  M  X2  —  w( 


Equating  the  coefficients  of  sinXt  and  t  on  each  side  gives  two  simultaneous 
equations  for  the  determination  of  P  and  A, 

P  P  X 


0  =  v0  - 


KP2^A  =  +  P 


mX  M  A2  —  w2  ’ 
1 


mX2  M  A2  —  w\ 


2  • 


After  eliminating  P  we  obtain 


212 


A5  =  v0 


[(m  +  M)A2  —  Mwf] 
A3K3m?M2(X2  —  w2)2' 


(2.2.6) 


(2.2.7) 
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Then  P  is  determined  from  the  first  of  the  simultaneous  equations  [6]  which 
gives 

_  mMA(A2  —  w\) 

0  (m  +  M)  A2  —  Mw\ 

Numerical  results  obtained  in  [6]  showed  that  this  method  is  fairly  good  in  some 

cases,  but  may  cause  a  large  error  in  other  cases.  This  indicated  that,  the  error 

in  the  solution  lies  not  in  the  mathematical  approximations  adopted  but  in  the 

physical  assumption  that  only  the  fundamental  mode  of  the  beam  affects  the 

impact  process.  Hence  some  criterion  must  be  found  which  will  give  conditions 

justifying  such  an  assumption. 

Among  the  other  approximation  methods  found  in  the  literature  is  the  energy 
.method  devised  by  Zener  and  Feshbach  [14],  and  applied  by  Lee  [6]  to  central 
impact  of  a  sphere  on  a  simply  supported  beam. 

2.2.3  Energy  method 

An  important  contribution  to  the  understanding  of  impact  characteristics  was 
made  by  Zener  and  Feshbach  [14]  in  their  considerations  of  energy  transfer. 
The  procedure  in  the  application  of  the  principle  of  conservation  of  energy  is 
to  assume  a  contact  force  variation  in  terms  of  which  the  energy  imparted 
to  the  beam  can  be  represented.  The  vibrational  energy  absorbed  by  each 
mode  is  composed  of  the  sum  of  kinetic  energy  and  strain  energy  of  bending 
at  a  particular  time.  Recall  that  the  deflection  of  the  beam  can  be  written  as 
w(x*,  t )  =  Cn(t)Wn(x*)  by  using  the  Green’s  function  method,  where  Cn(t) 
is  modal  coordinate  and  is  governed  by  the  following  differential  equation, 

Cn(t)  +w2nCn{t)  =  Wn(x*)f  (t),  n  —  1,  2,  3  •  •  • 
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It  is  easy  to  show  that  the  energy  for  the  nth  mode  during  impact  is  equal  to 

AE,  =  +  wlC2„(t)\ 

=  \w»(x')2 1  /r/(T)e”»"<ir|2,  (2.2.9) 

L  Jo 

where  the  T  is  the  time  when  contact  loses. 

Errors  resulting  from  ignorance  of  the  contact-force  variation  /(r)  can  be 
largely  eliminated  by  expressing  /(r)  in  the  form  of  a  normalized  force.  The 
rebound  velocity  of  the  impactor  can  be  expressed  in  the  form  ev0  by  introducing 
a  restitution  coefficient  e.  The  total  change  in  the  momentum  of  the  mass 
produced  by  the  impact  force  is  then 


fT 

mv0(  1  +  e)  =  /  f{r)dT. 
Jo 


(2.2.10) 


Now,  defining  the  normalized  impact  force  F(t)  by  the  relation  /(r)  =  mv0(  1  + 
e)F(r),  the  momentum  equation  becomes 

rT 


[  F(r)dT  =  1. 
Jo 


(2.2.11) 


Substitute  for  /(r),  equation  (2.2.9)  becomes 

A En  =  l-m2vl{  1  +  e)2Wn(x*)2\  F  F{r)eiw-Tdr\2.  (2.2.12) 

Li  J  0 

The  total  vibrational  energy  of  the  beam  is  equal  to  the  sum  of  the  energies  in  the 
separate  modes  of  vibrations,  i.e.,  the  total  vibration  energy  A E  =  A En. 

Let  us  define  a  function  R  as 


oo  fT 

R  =  mYj  Wn{x*)2\  /  F{T)eiWnTdr\2.  (2.2.13) 

n=l 

The  total  vibrational  energy  A E  now  reduces  to 

OO  1 

A E  =  ^2  A En  —  -mv q(1  +  e)2R  =  E(  1  +  e)2R.  (2.2.14) 

n=l  ^ 
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where  E  represents  the  initial  kinetic  energy  of  the  impactor.  Since  energy  is 
conserved,  AE  is  also  the  loss  of  energy  of  the  impactor,  giving  the  equation 

A E  =  E(  1  -  e2).  (2.2.15) 


Now  A  E/E  can  be  eliminated  from  equations  giving  the  coefficient  e  in  terms 
of  the  function  R ,  or 

e  =  ™  ■  (2.2.16) 


If  we  assume  that  an  expression  for  F(r)  satisfies  the  condition  (2.2.11),  then 
R  and  e  can  be  determined.  However,  the  F(t)  is  generally  unknown,  hence  an 
approximation  must  be  selected.  In  order  to  avoid  heavy  computational  burden 
to  evaluate  function  R,  a  simple  force  function  can  be  chosen.  One  example 
given  in  [6]  is  of  the  form 

F(t)  =  ^-sin(Trt/T0) 

=  -^-sinXt,  Vte[0,r0].  (2.2.17) 


where,  T0  is  the  estimated  impact  duration  by  assuming  two  rigid  body  impact, 
calculated  from  [5].  Hence,  the  impact  force  is  given  by 


/(t)  =  KosinXt. 


(2.2.18) 


where  the  K0  is  a  constant. 


2.3  Control  of  Impact  Forces 

As  we  have  pointed  out  in  chapter  1,  the  force  and  constrained  motion  control  in 
various  robotic  manipulations  is  the  main  motivation  to  this  research.  In  general, 
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there  are  three  modes  of  operation  so  that  a  robot  can  complete  a  task:  motion 
in  free  space,  impact  and  constrained  motion.  Obviously  the  manipulator  must 
change  from  one  mode  to  the  other  readily.  Usually  switching  from  the  con¬ 
strained  space  to  free  space  presents  no  problems.  However,  the  opposite  switch 
has  the  significant  problem  of  impact  forces.  The  impact  forces  can  be  very 
large,  and  can  drive  an  otherwise  stable  controller  into  instability.  Typically, 
it  is  the  force  control  strategy  that  must  deal  with  this  transient  phenomenon. 
However,  the  natural  elasticity  of  impacting  bodies  (demonstrated  by  the  nonlin¬ 
ear  impact  model),  can  cause  the  manipulator  to  rebound  from  the  environment. 
Thus,  the  manipulator  is  once  again  unconstrained.  This  phenomenon  can  cause 
oscillatory  behavior  [43].  Obviously  it  is  the  goal  of  any  controller  to  pass  this 
transitory  period  successfully,  and  have  the  manipulator  stably  exerting  forces 
on  the  environment  in  the  end.  The  controller  must,  therefore,  pass  through  the 
impact  phase  by  attempting  to  maintain  contact  with  the  environment  until  all 
of  the  energy  of  the  impact  has  been  absorbed.  Most  early  research  treated  the 
impact  as  a  transient  and  did  not  account  for  the  impact  effect,  which  resulted 
in  performance  tradeoffs.  This  problem  has  drawn  more  attention  recently  since 
high  precision  control  of  robotic  manipulators  has  become  increasingly  important 
in  a  variety  of  industrial  applications,  e.g.  laser  beam  technology,  semiconductor 
wafer  manufacturing  etc.,  thus  requiring  extra  attention  to  the  usual  dynamical 
effects  as  well  as  taking  into  consideration  otherwise  ignored  features  such  as 
dynamical  effects  due  to  impact.  Below  we  describe  some  the  impact  control 
schemes  that  have  been  developed  in  the  past  few  years. 
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2.3.1  Open-loop  control  strategy 

One  useful  method  that  has  been  used  for  controlling  impact  dynamics  is  open- 
loop  control  schemes.  The  motivation  for  this  strategy  stems  from  the  obser¬ 
vation  that  feedback  signals  may  obscure  effects  due  to  impact  because  of  the 
bandwidth  limitations  of  actuators  and  sensors  [4].  Salmatjidis  [3]  and  Chapnik 


Figure  2.2:  A  single-link  flexible  robot 

etc.  [4]  independently  studied  impact  and  force  control  problem  of  a  single-link 
flexible  robot  manipulator  shown  in  Figure  2.2,  where  a  flexible  robot  with  a 
flexible  beam  is  mounted  on  a  moving  base  driven  by  a  torque  T(f).  Suppose 
the  flexible  robot  is  transversely  collided  with  an  impactor  with  mass  nn  and 
velocity  u0.  If  we  just  consider  the  transverse  deflection  of  the  beam  and  ignore 
the  rotary  inertia  effects,  by  the  extended  Hamilton  principle [44],  the  dynamic 
equation  of  the  beam  is 

cP“u  dv 4 

P-Qj. f  +  EIfaA  =  o  <  x  <  l,  (2.3.1) 

where  y(x,t )  =  w(x,t)  +  x6(t),  w(x,t)  is  the  deflection  of  the  beam,  and  9(t) 
is  the  rotational  angel  of  the  hub,  and  f(x*,t )  is  the  distributed  load.  Equa¬ 
tion  (2.3.1)  describes  an  infinite  dimensional  system.  Typical  finite  dimensional 
model  approximations  are  used  for  the  controller  design. 
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Finite  Element  Model: 

The  elastic  deflection  w(x,t)  is  to  be  approximated  using  finite  element  tech¬ 
niques  [44].  Assembling  the  element  mass  and  stiffness  matrices  leads  to  the 
semi-discretized  equations  of  motion  for  the  flexible  arm  systems, 

MQ  +  CQ  +  KQ  =  F  (2.3.2) 

where  M,  C,  and  K  are  the  global  mass,  damping,  and  stiffness  matrices,  respec¬ 
tively,  F  is  the  external  force  vector  (applied  torque  is  specified  in  the  first  entry 
of  this  vector),  and  F  contains  the  unknown  reaction  (impact)  forces  at  each 
element.  Q  =  [9,  uq,  wi,  w2,  w2 . ,  w n,  wn]t  for  some  finite  integer  N  >  0. 

Simulation: 

In  order  to  calculate  the  time-domain  response  of  the  beam  from  the  equation 
(2.3.2),  the  Newmark-/?  integration  scheme  [45]  has  been  implemented,  and  yields 
the  vector  Q,  Q ,  and  Q  from  equation  after  every  integration  time  step.  In  order 
to  solve  the  equation  completely,  a  model  of  impact  force  must  be  given. 

Impact  Model  and  Control: 

Salmatjidis  used  a  stereomechanical  impact  model  by  assuming  that  the  impact 
duration  is  negligible.  The  impulsive  force  is  defined  as 

f=tki  i(t)dt  <2-3'3) 

Plugging  the  impulsive  force  model  into  the  dynamic  equation  (2.3.2)  and  taking 
the  limit  process,  the  problem  is  left  to  determine  to  the  amplitude  of  /.  Based 
on  this  impact  model  the  authors  realized  that  a  contact  problem  is  actually  a 
problem  with  holonomic  constraints.  The  primary  kinematic  axiom  of  a  contact 
problem  is,  that  the  structures  involved  do  not  penetrate  each  other.  Hence 
the  kinematic  constraint  provides  an  extra  equation  necessarily  to  determine  the 
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impact  force  /.  An  open-loop  control  algorithm  is  proposed,  which  operates  the 
three  modes  of  the  flexible  robot  as  follows; 

1) .  The  controller  monitors  contact  through  the  FSR  (  force  sensitive  resister  ) 
output,  while  in  free  space. 

2) .  When  contact  is  detected,  an  open-loop  control  scheme  is  implemented. 
The  desired  displacement  and  velocity  profiles  that  the  system  should  follow 
are  specified  by  using  the  simplified  dynamic  model  derived  through  kinematic 
(holonomic)  constraints  which  hold  when  contact  is  established,  then  the  re¬ 
quired  torque  is  computed. 

3) .  Once  the  transient  phenomenon  is  over  and  contact  is  established,  the  force 
regulator  is  turned  on  (  PD  control),  making  the  last  minor  corrections. 
Experimental  results  showed  that  this  control  scheme  is  better  than  convention 
force  control  methods  without  accounting  for  the  impact  forces.  As  pointed  out 
by  the  authors,  this  method  is  still  in  a  quite  primitive  state. 

Chapnik  etc.  [4]  proposed  another  open-loop  approach.  An  impact  model 
(energy  method)  developed  by  [6,  5]  is  used  for  their  analysis  of  impact  force, 
where  the  finite  duration  of  impact  is  accounted  for 

f(t)  =  K0sinXt.  (2.3.4) 

The  control  strategy  is  based  on  a  computed  torque  scheme  proposed  by  Bayo 
[46].  As  Bayo  has  shown,  only  one  degree  of  freedom  needs  be  specified  from 
equation  (2.3.2)  in  order  to  solve  for  the  required  torque.  Once  this  variable 
has  been  specified,  by  an  equation  in  time  domain,  the  other  degrees  of  freedom 
can  be  calculated.  It  was  decided  that,  in  accordance  with  Bayo,  the  variable 
specified  would  be  the  global  tip  acceleration,  jijv.  This  allows  the  specification  to 
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include  the  final  resting  place  of  the  beam,  and  to  force  minimum  tip  vibration  in 
the  global  sense.  Once  the  required  variable  is  specified,  Fast  Fourier  Transforms 
are  used  to  calculated  the  required  torque. 

The  weakness  of  the  open-loop  control  schemes  is  obvious.  Without  any 
feedback  the  performance  of  the  control  is  determined  by  the  accuracy  of  the 
system  and  impact  models.  Modeling  errors  may  result  in  even  greater  impact 
forces  than  without  control. 

2.3.2  Proportional  explicit  force  control 

The  another  successful  method  is  proposed  by  Volpe  [9,  10],  an  experimental 
approach  devised  to  solve  impact  control  problems.  By  using  an  explicit  force 
control  and  properly  tuning  the  feedback  gain,  the  author  showed  that  impact 
dynamics  can  be  compensated  satisfactorily.  To  better  understand  this  control 
scheme,  let  us  consider  a  generic  force-based  explicit  force  controller  shown  in 
Figure  2.3,  where  G  is  the  plant,  H  is  the  controller,  and  R  is  the  feedfor- 


Figure  2.3:  Explicit  force  control  block  diagram 

ward  transfer  function,  and  L  is  a  force  feedback  filter.  The  plant  G  may  be 
represented  by  the  fourth  order  model  or  the  reduced  second  order  model  [43]. 
The  controller  H  is  usually  some  subset  of  PID  control  (  e.g.  P,  PI,  PD).  In 
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[9],  authors  implemented  a  force-based  explicit  force  controller  with  proportional 
gain  and  extra  feedback  for  reaction  force  compensation.  The  plant  G  has  be 
expanded  into  its  components.  It  has  been  shown  that  the  system  is  equivalent 
to  one  shown  in  Figure  2.4.  An  analysis  using  the  root  locus  method  showed 


Figure  2.4:  Impact  control  block  diagram 

that  the  stability  of  the  closed  loop  system  is  guaranteed  for  H'  >  0,  i.e.  nega¬ 
tive  proportional  force  control  gains  down  to  —1  are  stable.  The  experimental 
results  showed  that  H  -»  —  1  is  desirable  for  impact  control.  An  interpretation 
in  [10]  goes  as  as  follows; 

By  using  this  control  scheme,  the  controller  does  not  utilize  the  force  error  sig¬ 
nal,  since  H'  =  H  +  1  «  0.  However,  the  reaction  force  of  the  impact  is  directly 
negated  by  a  feedback  signal.  Viewed  this  way,  the  impact  controller  does  not 
bounce,  because  the  oscillations  in  the  commanded  force  and  those  in  the  ex¬ 
perienced  force  are  equal  and  opposite.  Thus  the  surface  is  at  a  node  of  two 
interfering  pressure  waves.  No  net  force  means  no  net  acceleration.  Any  initial 
oscillation  is  damped  out  by  natural  and  active  damping.  The  drawback  of  this 
control  algorithm  is  that  the  best  feedback  gain  for  impact  dynamics  is  very  poor 
for  tracking  purposes  which  results  in  performance  tradeoff. 
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2.3.3  Other  approaches 

There  are  some  other  approaches  aiming  to  control  impact  forces. 

A)  Passive  Compliance  and  Damping 

One  proposed  method  of  dealing  with  the  impact  problem  is  to  use  passive  com¬ 
pliance,  either  on  the  end  effector  or  in  the  environment.  Some  researchers  have 
proposed  the  use  of  soft  force  sensors  [47].  These  methods  appear  to  provide 
stable  impact  in  two  ways.  First,  the  material  used  naturally  provides  passive 
damping  that  helps  to  absorb  some  of  the  energy  of  impact.  Second,  the  compli¬ 
ance  of  the  material  effectively  lessens  the  stiffness  of  the  system  composed  of  the 
material  and  the  environment.  There  are  problems  with  the  passive  compliance: 
it  may  not  be  modified  without  physical  replacement  of  the  material  and  it  also 
limits  the  effective  stiffness  of  the  manipulator  during  position  control  etc. 

B)  Active  Damping 

Another  method  is  to  employ  maximal  damping  during  the  impact  phase  [48]. 
The  goal  of  this  strategy  is  to  damp  out  the  oscillations  caused  by  the  tran¬ 
sition.  While  this  may  be  successful  for  soft  environments,  stiff  environment 
may  present  some  problems,  since  the  feedback  signal  seems  obscure  the  effect 
of  impact  forces  [4]. 
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Chapter  3 


MODELING  OF  IMPACT  ON  FLEXIBLE 
STRUCTURES 


In  this  chapter,  we  will  provide  a  systematic  way  to  model  the  impact  forces  on 
flexible  structures  via  use  of  Hertz  law  of  impact.  For  the  sake  of  simplicity,  we 
only  consider  a  flexible  structure  subject  to  impact  forces  occurring  from  contact 
with  an  elastic  body.  Here  a  flexible  structure  means  both  finite-dimensional  and 
infinite-dimensional  systems  (  a  lumped-parameter  structure  and  a  distributed- 
parameter  structure) .  We  restrict  attention  to  the  problem  of  modeling  impact 
dynamics,  existence  of  solutions  to  the  model,  and  application  examples.  Nu¬ 
merical  aspects  will  be  addressed  in  chapter  3. 

In  section  2  of  this  chapter  an  impact  is  formulated  and  the  Hertz  equation 
is  derived  through  the  Hertz  law  of  impact.  In  section  3  we  will  discuss  some 
basic  properties  of  the  Green’s  function  associated  with  the  Euler  Bernoulli  beam 
equation.  This  equation  is  used  to  describe  the  motion  of  the  beam.  In  section 
4  we  establish  the  existence  and  uniqueness  of  solutions  of  the  Hertz  equation 
by  applying  the  contraction  mapping  principle.  We  will  discuss  two  examples  in 
section  5. 
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3.1  Formulation  of  the  Impact  Problem 


3.1.1  Impact  on  a  lumped-parameter  structure 

We  first  study  the  impact  between  an  elastic  body  and  a  lumped-parameter 
structure.  For  our  purpose,  the  impact  problem  can  be  formulated  as  follows; 
A  lumped-parameter  structure  is  described  in  Figure  2.1,  where  is  the  mass 
of  the  ith  rigid  body  and  ki  represents  the  stiffness  ( i  =  1, 2,  •  •  • ,  n  —  1,  n)  of  a 
spring  between  bodies  i  and  i  +  1.  The  words  ’’rigid  body”  certainly  need  to 
be  clarified  here,  the  meaning  of  ’’rigid  body”  is  macroscopic.  For  the  impact 
problem,  ’’rigid  body”  is  still  assumed  to  be  locally  deformable.  Suppose  that 
the  system  is  struck  by  a  mass  m  having  a  spherical  surface  at  the  point  of 
contact  with  initial  moving  velocity  v0.  We  further  assume  that  the  Hertz  law 
of  impact  is  valid  ,  i.e, 

a  =  (3.1.1) 

where,  a(t),  the  relative  approach,  is  the  difference  between  the  displacements 
of  the  first  rigid  body  and  the  contacting  body,  measured  from  the  instant  of 
initial  contact.  Hence, 


a  =  s(t )  —  xx  (f),  (3.1.2) 

where  x\  (t)  is  the  displacement  of  the  first  rigid  body  of  the  system,  s(t)  is  the 
displacement  of  the  ball  under  the  contact  force  /(f),  and  is  given  by, 

s(t)  =  v0t  -  —  f  /(r)(t  -  r)dr.  (3.1.3) 

m  Jo 
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Figure  3.1:  Impact  on  a  lumped-parameter  structure 

From  equations  (3.1.1),  (3.1.2)  and  (3.1.3),  we  obtain  the  following  nonlinear 
integral  equation, 

K[f(t)}2/ 3  =  v0t  -  —  [  f(r)(t  -  T)dT-xi(t).  (3.1.4) 

m  Jo 

If  we  assume  that  the  flexible  structure  is  at  rest  just  before  impact,  it  is  easy 
to  show  that  the  displacement  x\(t)  can  be  expressed  as 

xi(t)  =  [  Gi(t  -  T)f(r)dT  (3.1.5) 

where  the  function  Gi(-)  is  the  Green’s  function  of  the  system.  Replacing  x\ (t) 
in  equation  (3.1.4)  by  (3.1.5),  we  get  the  nonlinear  integral  equation 

K[f(t)]2/3  =  vot  -  —  [  f{r){t  -  r)dr  —  [  Gi(t  -  r)/(r)dr.  (3.1.6) 

m  Jo  Jo 

We  note  here  that  in  general  it  is  impossible  to  solve  (3.1.6)  analytically. 
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3.1.2  Impact  on  a  flexible  beam 


We  now  study  the  impact  between  an  elastic  body  and  a  flexible  beam  with 
various  boundary  conditions.  Suppose  a  beam  is  struck  transversely  by  a  mass 
m  having  a  spherical  surface  at  the  point  of  contact  with  initial  moving  velocity 
wo  (  see  Figure  2.2).  In  the  same  way  as  for  the  lumped-parameter  structure  in 


Figure  3.2:  Impact  on  a  flexible  beam 

section  2.2.1,  we  end  up  with  the  nonlinear  integral  equation 

K[f{t)]2/ 3  =  v0t  -  —  [  f{r)(t  -  r)dr  -  w(x*,t).  (3.1.7) 

m  Jo 

where  w(x*,t )  is  the  deflection  of  the  beam  at  the  point  of  contact  x*.  Before 
we  can  carry  out  further  analysis,  it  is  necessary  to  represent  the  deflection  of 
the  beam  at  the  point  of  contact  w(x*,t). 

3.2  Deflection  of  the  Beam  and  Green’s 
Function 

If  we  restrict  attention  to  transverse  vibrations  only  and  assume  that  the  beam 
is  long  and  slender,  the  transverse  shear  and  torsional  effects  can  be  neglected, 
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and  the  dynamics  of  the  beam  can  be  described  by  the  Euler  Bernoulli  beam 
equation.  The  deflection  of  the  beam  is  in  turn  obtained  by  solving  the  following 
partial  differential  equation 

d^w  - 

p-^w  +  El  A  w  =  f(x,  t)  0  <x<l,  (3.2.1) 

otz 

in  which  A  =  l  is  the  length  of  the  beam,  p  is  the  mass  density  and  El 
is  the  bending  stiffness  (here  p  and  El  are  assumed  constants),  and  f{x,t)  is 
the  distributed  load.  The  deflection  of  the  beam  is  uniquely  determined  when 
equation  (3.2.1)  is  solved  under  the  appropriate  initial  and  boundary  conditions. 
If  we  assume  that  the  beam  is  at  rest  just  before  impact,  the  initial  conditions 
are 


w(x,0)  =  w(x,  0)  =  0.  (3.2.2) 

Various  boundary  conditions  of  interest  can  be  described  as 

Biw(x,t)  =  0,  *  =  1,2,  x  =  0,1,  (3.2.3) 

where  Bi  is  a  linear  homogeneous  differential  operator  of  maximum  order  3. 

The  concentrated  load  f(t)  is  obtained  as  a  limiting  process  of  a  uniformly 
distributed  load  f(x,t )  over  a  small  range  25  of  the  beam.  Thus,  by  letting 
f(x,  t)  —¥  oo  while  5  — »  0,  the  contact  force  f(t)  is  obtained  as 

f(t)  =  lim  [  f(x,t)dx  (3.2.4) 

5— >0,f— >oo  Jx‘—  <5 

One  of  the  most  popular  methods  for  analyzing  linear  partial  differential  equa¬ 
tions  is  the  Green’s  function  method  [49].  With  the  aid  of  a  Green’s  function, 
the  solution  of  a  certain  class  of  PDE  can  be  expressed  as  an  integral. 
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Definition  3.2.1  A  function  G(x,(',t)  €  L2[0,l]  is  called  a  Green’s  function  of 
(3.2.1)-(S.2.3),  if  it  satisfies  the  following  conditions: 

i)  As  a  function  of  the  argument  x,  it  satisfies  the  homogeneous  differential  equa¬ 
tion  ,  i.e.  f  —  0,  everywhere  except  at  x  —  (  where  it  may  have  a  singularity. 

ii)  As  a  function  of  the  argument  t,  it  satisfies  the  homogeneous  differential  equa¬ 
tion  everywhere  except  at  t  —  0  where  it  may  have  a  singularity. 

Hi)  As  a  function  of  the  argument  x  ,  G(x,(]t)  satisfies  the  boundary  condition 
(3.2.3). 

iv)  It  satisfies  the  initial  conditions 

G{x,C,  0+)  =  0  and  dG(xff),0  )  _  ^x  (-yp  (3.2.5) 

Note  that  (i)  and  (ii)  above  lead  to 

{EI  A  +p~}G(x,C,t)  =  Hx.om-  (3.2.6) 

If  the  PDE  (3.2.1)-(3.2.3)  is  self-adjoint  (which  depends  on  the  nature  of  the 
boundary  conditions).  Cantilevered  and  simply  supported  beams  result  in  self- 
adjoint  PDE’s.  G(x,ff,t)  is  symmetric  with  respect  to  x  and  (,  and  can  be 
expressed  in  terms  of  an  eigenfunction  expansion.  Therefore,  in  this  case  it  can 
be  shown  that  the  Green’s  function  for  the  PDE  (3.2.1)-(3.2.3)  can  be  expressed 
in  the  form  (we  refer  the  reader  to  [50]  for  details) 

G(x,G,t)  =  f ]Wk(x)Wk(0S-^^H(t),  (3.2.7) 

*=i  w * 

where  H{t)  is  the  unit  step  function,  {Wk(x)}'f(Ll  is  an  orthonormal  basis  of 
eigenfunctions  and  {u;fc}^T1  are  the  corresponding  eigenvalues.  The  following 
theorem  is  standard  in  the  theory  of  partial  differential  equations. 
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Theorem  3.2.1  [50]  A  representation  of  the  solution  to  the  PDE  (3.2.1)-(3.2.3) 
in  terms  of  the  Green’s  function  G(x,  t)  is 

w(x,t)  =  [  [  G(x,G,t-T)f(C,r)d([dT  (3.2.8) 

Jo  Jo 

□ 

For  the  impact  problem,  because  the  contact  can  be  treated  as  a  point  contact, 
the  contact  force  has  the  special  form  (3.2.4).  Hence  equation  (3.2.8)  can  be 
further  simplified  as 

w(x,t)  =  f  G(x,x*]t  —  r)/(r)dr.  (3.2.9) 

Jo 

For  simplicity,  we  write  G(x*]t)  instead  of  G(x*,x*;t)  in  the  rest  of  the  disser¬ 
tation.  From  equations  (3.1.7)  and  (3.2.9),  the  Hertz  equation  becomes 

K[f(t)}2/3  =  v0t-—f  f(r)(t-T)dr  -  f  G(x*;t  -  r)f(r)dr.  (3.2.10) 
m  Jo  Jo 

Remark  3.2.1  For  the  impact  problems  in  section  3.2.1  and  3.2.2,  we  have 
derived  the  equations  of  impact  dynamics  (3.1.6)  and  (3.2.10)  by  means  of  the 
Hertz  law  of  impact.  Note  that  both  equations  have  the  same  structure;  the  only 
difference  is  in  the  nature  of  their  Green’s  functions.  For  a  finite- dimensional 
system,  the  Green’s  function  is  much  simpler  in  general,  we  only  use  equation 
(3.2.10)  for  analysis. 

3.3  Analysis  of  the  Hertz  Equation 

We  should  point  out  that  there  already  exist  some  equations  similar  to  equation 
(3.2.10)  to  model  impact  dynamics  on  flexible  structures  [5,  6].  Timoshenko  [11] 
derived  one  for  a  simply  supported  beam.  Though  these  equations  haven’t  been 


34 


analyzed  in  any  great  detail  in  the  existing  literature,  some  numerical  methods 
for  solving  these  equations  have  been  presented  in  some  detail.  Our  view  is 
that  some  theoretical  analysis  is  necessary  for  both  proving  the  validity  of  this 
equation  and  developing  efficient  numerical  methods.  T  contraction  mapping 
technique  is  employed  here  to  show  that  a  unique  solution  exists  for  the  Hertz 
equation.  Before  invoking  the  contraction  mapping  theorem  some  rearrange¬ 
ments  are  necessary.  Let 

L(t )  =  t  +  mG(x*]t).  Vi  >  0  (3.3.1) 

Equation  (3.2.10)  can  be  rewritten  as 

f{t)2'3  =  v'0t - —  f  f(r)L(t  -  r)dr,  (3.3.2) 

m  Jo 

where  v'Q  =  v0/K]  m!  —  mK, 

f(t)  =  Wot  -  “  Jo  f(r)L(t  -  r)dr]3/2  (3.3.3) 

=  K*  “7  L  fW)(t  -  T)dT~\?  f  fW)G(x *;  t  -  r)dr]3/2.  (3.3.4) 
m  Jo  K  Jo 

Note  that  v'0  is  assumed  to  be  positive  always.  Both  equations  (3.3.3)  and  (3.3.4) 
will  be  used  in  the  following  analysis.  The  following  result  is  well  known  (see, 
e.g.,  [51]). 

Theorem  3.3.1  (Contraction  Mapping  Theorem)  Let  X  be  a  Banach  space,  and 
B  a  closed  subset  of  X.  Let  P:  B  — »  B  be  an  operator  satisfying  the  following 
condition:  3  p  <  1  such  that 

\\Px-Py\\  <  p\\x  —  y\\,  yx,y£B. 

Then 

a)  P  has  exactly  one  fixed  point  in  B  (denoted  by  x*). 
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b )  For  any  x0  €  B,  the  sequence  {xn}o°  defined  by 


xn+i  -  Pxn,  n  >  0 

converges  to  x*.  Moreover, 

Pn 

*»-**  <  z - \\Px0  -  loll- 

1  -  p 

We  will  use  this  theorem  as  the  main  tool  to  show  that  equation  (3.3.3)  has 
a  unique  solution  by  constructing  a  contraction  operator  P  on  an  appropriate 
closed  subset  B  of  a  Banach  space. 

3.3.1  Local  existence  and  uniqueness 

We  shall  first  establish  some  conditions  under  which  (3.3.3)  has  exactly  one  so¬ 
lution  over  every  finite  interval  [0,5]  for  6  sufficiently  small,  i.e.,  conditions  for 
local  existence  and  uniqueness.  We  shall  then  obtain  stronger  conditions  for 
global  existence  and  uniqueness,  i.e.,  conditions  under  which  (3.3.3)  has  exactly 
one  solution  over  [0,  T]  for  some  finite  T.  Hypotheses  in  both  theorems  are  sat¬ 
isfied  by  finite  dimensional  systems  as  well  as  simply  supported  and  cantilevered 
beams. 

Theorem  3.3.2  Suppose  that  the  Green’s  function  G(x*;  t)  is  uniformly  bounded 
over  x*  €  [0,  l ]  and  t  6  [0,  T]  for  some  T  >  0.  Then  there  exists  a  small  enough 
8  >  0  such  that  (3.3.3)  has  a  unique  continuous  solution  for  t  e  [0,5]. 

Proof  :  Let  M  >  0  be  such  that 

\G{x*;t)\  <  M  Vt>  0  and  Vx*e[0,i]. 

Let  N  >  0  be  a  sufficiently  large  constant.  Let  5  >  0  be  small  enough  such  that 

M  4  < 
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fa)  $  <  _ _ • 

[UJ  °  ^  [l/m'+M/K]  ’ 

(in)  2 (52/m'  +  M8/K)yJ{v'Q  +  MN/K)6  <  1. 

Lei  i/ie  Banach  space  be  C[0,  (5],  the  space  of  continuous  real  valued  functions 
from  [0,5],  endowed  with  the  sup  norm,  i.e.  ||/||oo  =  supiej0)(5]  |/(i)|.  Let  us 
define  the  mapping  P  :  C[0, 5]  — >  C[0, 5]  by 

1  rt  3/2 

Pf{t)  =  Wot  -  —  j(  f(T)L(t  -  r)dr ]  ,  Vi  €[0,5].  (3.3.5) 

The  domain  of  P  is  defined  by 


BM=  {/(•)€  C[0,  i];JV>  /(f)  >0; 

vot  ~  J  f(T)L(t  -  r)dr  >  0  Vi  G  [0,5]}.  (3.3.6) 

Obviously,  Lf[0, 5]  is  a  closed  subset  of  the  Banach  space  of  continuous  functions 
on  [0, 5]. 

The  rest  of  the  proof  is  divided  into  two  parts:  first,  we  show  that  P  maps 
L?[0,5]  into  itself.  Then  we  show  that  P  is  a  contraction  mapping  on  B[0, 5]. 
a)  Pf  >  0  V/  G  B[0, 5]  by  definition.  Define  mapping  F  :  B[0, 5} B[0, 6} 
by 


Ff(t)  =  K*  ~  J0  f(T)L(t  ~  T)dr] 

Now  f  G  5[0, 5]  implies  ^  fo  f(r)(t  —  r)dr  >  0.  Hence, 

\Ff(t)\  =  \Wot~  Jo  f(T)(t~T)dT  -  ^  Jo  f(r)G(x*;t-T)dT]\ 

<  v'ot  +  ^JQ  \f(r)\\G(x*',t  ~  r)\dr 

<  (v'0  +  MN/K)t,  t  G  [0,5], 


Ff(t)  >  0  by  definition.  Therefore, 


F f(t)  <  (v'  +  MN/K)t 

implies  Pf{t)  <  [(^o  +  MN/K)t]3/2  <t,  Vi  G  [0,5], 
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FPf(t) 


> 

> 


V'ot~m/L  Pf(TW-T)dT~Jt  J0  Pf(r)G(x*]t-r)dr 

1  rl  1  /•« 

v'0t - -  /  (t  —  t)t<It  —  —  rMdr 

m  Jo  K  Jo 

,  1  t 2  t2M 

Vot  ~  m'2  ~  ~2K 


>  0  We  [0,5]. 


Thus,  we  have  shown  that  PB[(),  6]  c  B[0,  5]. 

b)  V/,,/,€B[0,i|, 

Pfi(t)  -  Pfi{t)  =  Wot  Jo  h(r)L(t  -  r)rfr]3/2 

-K*  -  ^7  JQ  f2(r)L(t  -  r)dr]3/2 

Let  x(t)  =  Wot  ~  ^  Jo  fi(T)L{t  -  r)dr]1/2. 

y(t)  =  Wot  Jq  f2(r)G(t  -  T)drW2. 

Since  fi,fs  €  -B[0, <5]  imply  x(t),y(t)  are  well  defined. 


Pfi(t )  -  PfsO) 
\Ph(t)~Pf2(t)\ 
x2(t)-y2(t) 


At)-y\t)  I 


=  x3(t)-y3(t)  Vie  [0,5] 

<  \x2(t)-y2(t)\\x(t)+y(t)\ 

=  vot  ~~j  h  (r)L(i  -  r)dT 

-Wot  jQ  h (T)L(t  -  T)dr) 

=  ~i  [  (MW  -  fi(r))L(t  -  r)dT. 

fib  J  L) 

=  1^7^  U2W)  -  fi{r))L{t  -  r)dr\ 

<  ^  Jo  l/2(r) -/i(r)|(i-r)dr 

+1<  Jo  l^2(r)  ~  fiW)\\G{x*-,t  -  T)\dT 

<  W/m'  +  Mt/K)  H/2-/1II00 

<  (62/m'  +  M6/K)\\f2-f1\\00. 
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\x(t)+y(t)\  <  \x(t)\  +  \y(t)\  <2y/v'0t  + MNt/K 

<  2^/  (vq  +  MN/K)5. 

\Pfi(t)-Pf2(t)\  <  \x2(t)-y2(t)\\x(t)  +  y(t)\ 

<  2(62/m'  +  M6/K)y/(v'0  +  MN/K)6\\f2(-)  -  /i(-)IU 

<  p||/2  —  /i||oo  We  [0,6]. 

where  p  —  2(<52/m'  +  MS/K)yJ(v'0  +  MN/K)5,  and  by  the  property(iii)  of  8, 
p  <  1.  Therefore, 

HP/i-JValleo  =  sup  \Pfi(t)  -  Pf2(t)\ 

te[o,8] 

<  p||/2  —  /llloo- 

so  that  P  is  a  contraction  mapping  on  5[0,  <5]. 

Finally,  using  theorem  2-4- 1,  it  follows  that  the  mapping  P  has  a  unique  fixed 
point  in  5[0,<5].  It  is  clear  that  f  is  a  solution  of  the  equation (3. S. 3)  over  [0,(5] 
if  and  only  if  Pf  =  f ,  i.e.  f  is  a  fixed  point  of  P  over  J3[0,  <5].  This  completes 
our  proof.  □ 

The  above  theorem  shows  that  a  unique  solution  exists  over  t  €  [0, 5]  for  some 
small  (5.  Our  interest  is  to  find  the  impact  force  variation  during  the  entire 
contact  period.  The  following  theorem  will  establish  this  global  result. 

3.3.2  Global  existence  and  uniqueness 

Theorem  3.3.3  Suppose  that  equation  (3.3.3)  has  a  local  unique  solution  over 
[0, 5]  for  some  sufficiently  small  6.  If  f(5)  >  0,  then  3e  >  0,  such  that 
equation(3.3.3)  has  a  unique  continuous  solution  on  [0,5  +  e]. 
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Proof  ;  The  argument  is  similar  to  the  proof  of  the  local  version.  We  will  only 
carry  out  the  details  of  a  crucial  step  here. 


Let  g  :  [0,  <5]  —>  R  be  the  unique  solution  of  (3.3.3)  established  in  the  proof 
of  theorem  2. f.2.  Let  N  be  a  positive  number  larger  than  ||<?||oo-  Let  e  >  0  be  a 
small  positive  constant.  Let 

D[ 0,  (5  H-  e]  =  {/  €  C[0,  <5  +  e];  /  |[o,a]=  g\ 

v'0t  — [  /(r)L(i  -  r)dr  >0,  Vi  £  [0, 6  +  e]}.  (3.3.7) 
m  Jo 

Clearly  D  is  a  closed  subset  of  (C[0,  S  +  e],  ||  •  ||oo). 

Let  F  \  D[ 0, 6  4-  c]  — ^  C[ 0,  S  +  c]  be 

Ff(t)  =  Vo*  ~  fQ  /(rW  -  T)dr> 

and,  let  P  =  F3/2. 

We  will  show  that  for  small  enough  e,  P(D )  C  D,  and  P  is  a  contraction 
mapping,  thus  establishing  the  theorem.  Note  that  it  follows  easily  as  in  the  proof 
of  theorem  2.4.2  that 


\Pf(t)\  <  N  Vte[0,5  +  c\. 

The  crucial  step  is  to  show  that  FPf(t )  >0  Vi  6  [0, 5  +  e],  V/  6  D.  Now, 
FPf(t)  —  f(t)  Vi  <  5  since  f  |[0,$]  satisfies  (3.3.3).  For  0  <  i  <  e, 

FPf(t  +  6)  =  v'06-~  f  P f(r){t  +  5  -  r)dr 

m'  Jo 

So  Pf(T)G(x*’t  +  5  ~T)dT 

+v'ot  jQ  Pf(T  +  5W  ~  T)dT 

=  v'08  -  -1  [S  Pf(r)(5  -  r)dr  -  ±  [*  Pf(r)G(x *;  5  -  r)dr 

m  Jo  K  Jo 
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~1<  Jo  t  +  S~T)-  G(x*]5  -  r))dr 

^  ~bC  nf)dT  ~bllP!{T+ s)L{t  - T)dT 

=  [/2/3W  -  £7  Jo  f{r)dT  -  G(x*;t)  +  G(x*;0) 

+v'ot  ~  ^7 1  pfij  +  8)L{t  -  r)dr } 

where  G(x*; t)  =  /05  f(r)G(x*;t  +  5  —  r)dr.  By  the  continuity  of  the  Green’s 

function,  for  small  enough  e, 

\G(x-,t) -G(x\ 0)1  <  \f2/3(sy,  Vt  e  [0, e] 

For  such  e, 

FPf{t  +  S)  >  0  Vte[0,t]. 

Hence,  we  have  shown  that  PD[ 0,  <5  +  e]  C  -D[0, 5  +  e]. 

It  is  easy  to  show  that  P  is  a  contraction  mapping.  The  existence  of  a  unique 
solution  on  [0, 5  +  e]  follows  at  once  now.  □ 

Remark  3.3.1  The  global  version  has  the  following  physical  interpretation.  The 
condition  f(t)  >0  means  that  the  contact  is  in  progress  att  >  0;  finally,  f(T)  — 
0  means  that  the  objects  are  just  about  to  cease  to  be  in  contact,  i.e.,  T  is  the 
impact  duration. 

3.4  Two  Applications 

3.4.1  A  single-link  flexible  robot 

High  precision  control  of  robotic  manipulators  is  becoming  increasingly  impor¬ 
tant  in  a  variety  of  applications,  e.g.,  laser  beam  technology,  semiconductor  wafer 
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manufacturing  etc.  This  requires  paying  extra  attention  to  the  usual  dynamical 
effects  as  well  as  taking  into  consideration  otherwise  ignored  features  such  as  dy¬ 
namical  effects  due  to  impact.  Here  we  consider  a  flexible  robot  with  a  flexible 
beam  mounted  on  a  moving  base  driven  by  a  torque  T(t).  Suppose  the  flexible 
robot  transversely  collides  with  an  impactor  with  mass  m  and  velocity  v0.  If  we 
just  consider  the  transverse  deflection  of  the  beam  and  ignore  the  rotary  inertia 
effects,  by  the  extended  Hamilton  principle  [44],  the  dynamic  equation  of  the 
beam  is 

d2v  dv * 

+  EI~^  =  /(**»*)>  0  <  x  <  l.  (3.4.1) 

where  y(x,  t)  =  w(x,t)  +  x 9(t),  w(x,t )  is  the  deflection  of  the  beam,  and  6{t) 


Figure  3.3:  A  single-link  flexible  robot 

is  the  rotation  angel  of  the  hub  and  the  f(x*,t)  is  the  distributed  load.  For  the 
impact  problem,  the  impact  force  has  the  general  form 

f(t)  =  lim  f  f(x,t)dx  (3.4.2) 

6— >0,f— >oo  Jx*—5 

where  x*  is  the  contact  point  of  the  beam.  The  boundary  conditions  are 
dy 2  dv 3 

EIdx*  ~lHWdxt +T^  =  °’  ^)=0’  at  x  =  0  (3A3) 

dy2  dy 3 

EI-£-  =  El^r  =  0,  at  x  =  l. 
dxz  dxA 
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where  Ig  is  the  inertia  of  the  base  and  T(t)  is  the  torque  applied  to  the  hub. 

In  order  to  apply  the  Green’s  method  as  described  in  section  3.3  to  obtain 
the  deflection  of  the  flexible  robot,  we  first  need  to  solve  the  eigenvalue  problem 
of  equation  (3.4.1)-(3.4.3).  The  general  methods  can  be  found  in  [52].  Posbergh 
[53]  applied  spectral  operator  theory  to  carry  out  some  rigorous  analysis.  For 
the  sake  of  brevity  some  useful  results  are  listed  here  (  see  [53]  for  details) . 

The  dynamic  equations  are  first  put  into  the  infinite-dimensional  state-space 
form 

d 

— s(t)  =  As(t)  +  Bu(t), 
at 

where  s(t)  are  state  variables,  u(t)  are  inputs.  A  and  B  are  two  abstract  opera¬ 
tors  defined  on  some  Hilbert  space  H. 

Rl:  The  operator  A  is  a  skew  adjoint  operator.  ^ 

R2:  The  spectrum  of  operator  A  is  discrete  and  the  characteristic  equation  is 
I  3s 

— — (1  +  cos(/3l)cosh(/3l))  —  sin(fll)cosh(f3l )  +  sinh((3l)cos((3l )  =  0  (3.4.4) 
where  the  eigenvalues  A  are  given  by  A2  =  —  -y-/?4. 

In  general,  there  are  a  countably  infinite  number  of  f3  which  satisfy  this  equation, 
denoted  by  (3n,  n=  1, 2,  •  •  •,  giving  rise  to  a  discrete  set  of  eigenvalues. 

Remark  3.4.1  If  IH  =  0  (the  hinged-free  case),  the  equation ( 3. f. 4)  reduces  to 

sin(/3nl)cosh(/3nl)  =  sinh(/3nl)cos((3nl). 

Remark  3.4.2  If  Ih  — >  oo  (the  cantilevered  case),  the  equation  ( 3.  f.f)  reduces 
to 


cos((3nl)cosh(/3nl )  =  —  1. 
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R3:  The  associated  eigenfunctions  are  orthonormal  and  have  the  following  prop¬ 
erties 


K"(x) 

0  <  X  <  l, 

(3.4.5) 

K(  o) 

=  ~<K(  o), 

X  —  0. 

<Wn,Wm> 

—  dn,m 

where  the  inner  product  defined  on  H  is  given  by 

<  f,g  >=  [  pfgdx  +  IH 
Jo 


df(0)ig(0) 
dx  dx 


Next,  we  derive  the  Green’s  function  for  equations  (3.4.1)-(3.4.3).  Since  our 
interest  is  to  model  the  impact  dynamics,  we  will  set  T(t)  =0  for  all  t. 

Because  {W^}^  are  orthonormal  and  form  a  complete  basis  of  the  Hilbert 
space  H,  the  solution  of  equation  (3.4.1)  can  be  expressed  in  the  form  [44] 


y(x,t)  =  £  Wk(x)Ck(t), 


(3.4.6) 


fc=0 


where  Ck{t)  is  the  kth  modal  coordinate.  From  equations  (3.4.1)  and  (3.4.5), 


Ck(t)  +  w2kCk(t)  =  Wk(l)f(t ), 

{Ib  +  IhW)  =  0.  (3.4.7) 


If  we  assume  that  the  flexible  robot  is  at  rest  before  impact,  the  corresponding 
initial  conditions  are  zero,  and  the  solution  y(x,t )  can  be  expressed  as 

y(x,t)  =  [  G(x,x*-t  -  r)f(r)dT  (3.4.8) 

Jo 

where  the  Green’s  function  G  is  given  by 

k= 1  Wk 
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In  order  to  apply  theorems  2.4.2  and  2.4.3,  we  then  need  to  show  that  the 
Green’s  function  G(x,(]t)  is  uniformly  bounded.  Observe  that  from  R2  we  get 
two  boundary  cases  (hinged-free  and  cantilevered-beam  case)  as  the  inertia  of 
the  base  Ih  approaches  two  extreme  cases  (0  and  oo),  it  is  shown  in  Appendix 
A. 2  that  the  Green’s  functions  are  uniformly  bounded  in  the  two  cases.  It  is 
possible  to  carry  out  similar  analysis  for  the  other  cases,  but  we  will  not  do  so 
here. 

3.4.2  A  smart  actuator 

Piezoelectric  and  magnetostrictive  actuators  are  gaining  increasing  attention  for 
their  potential  use  as  positioners,  motors,  and  vibration  suppressors  [54,  55]. 
They  are  anticipated  to  play  a  prominent  role  in  high  precision  manufacturing 
devices  for  optical  instruments  such  as  lasers  and  cameras,  and  high  precision  po¬ 
sitioning  in  semiconductor  chip  manufacturing  chip  etc.  They  are  being  applied 
already  in  high  precision  optical  telescopes,  cameras  etc.  [56,  57]. 

An  actuator  using  piezoelectric  and  magnetostrictive  materials  is  being  de¬ 
signed  and  fabricated  by  R.  Venkataraman  [58]  as  a  part  of  a  project  on  smart 
structures  technology  with  applications  to  rotorcraft  systems  at  the  University 
of  Maryland.  A  description  of  the  actuator  is  shown  in  Figure  2.4,  where  a 
piezoelectric  stack  and  two  magnetostrictive  rods  are  electrically  connected  and 
controlled  by  a  single  sinusoidal  voltage  signal  v(t).  During  the  electrical  half 
cycle  the  piezoelectric  stack  expands  and  clamps  onto  a  disk,  while  the  magne¬ 
tostrictive  rod  pushes  the  clamp.  The  combined  motion  results  in  the  rotation 
of  the  disk.  During  the  second  half  cycle  the  piezoelectric  stack  releases  its  hold 
on  the  disk  and  the  magnetostrictive  actuators  return  to  their  initial  positions. 
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Figure  3.4:  A  sketch  of  smart  actuator 

The  speed  rotation  of  the  disk  is  regulated  by  the  frequency  of  the  power  supply. 
Phase  angle  relationships  of  the  piezoelectric  stack  and  magnetostrictive  rods  are 
given  in  Figure  2.5.  During  the  expansion  of  the  piezoelectric  stack,  it  collides 
with  the  disk  and  an  impact  occurs.  Proper  modeling  of  such  impact  dynamics 
is  important  because: 

1) .  The  impact  force  may  become  so  large  that  it  will  shatter  the  materials. 
The  impact  model  described  earlier  in  the  chapter  can  be  used  to  select  materi¬ 
als  (clamp  and  disk)  as  well  as  to  determine  a  proper  impact  velocity. 

2) .  The  impact  duration  is  non-zero  and  finite,  and  this  fact  turns  out  to  be  very 
important  for  design  considerations.  In  simulations  it  has  been  found  that  the 
average  impact  duration  could  be  as  high  as  50  times  the  time  step  in  adaptive- 
step  size  numerical  simulations  which  ignore  this  phenomena. 

The  disk  and  the  clamp  system  is  modeled  as  a  spring-mass  system  as  shown  in 
Figure  3.6.  We  assume  that  the  Hertz  law  of  impact  is  valid, 
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Figure  3.5:  The  relationship  of  expansion  and  contraction 


Force  due 
to  input  voltage 


Figure  3.6:  Disk  and  clamp  system. 


47 


K  [f  (()]2/s  =s(t)-w  (() , 


(3.4.9) 


where  w  ( t )  is  the  displacement  of  the  disk  at  the  point  of  contact,  s  (t)  is  the 
displacement  of  the  clamp  under  the  contact  force  /  (t),  and  is  given  by, 


^  ^  _  V  clamp  SlTl  {p)n  clampt)  C\V  StTl  (uJn  clampt) 

OJ n  clamp  ^  dampen  clamp 

1  /■< 

-  /  /  (r)  sin  [w  (t  ~  7")]  dr 

o  nirlamn  JO 


^ n  clamp  clamp 


where 


Ci  = 


de  Ae  Ne 
se  L 


and  w  ( t )  is  given  by 

/,\  Vdisk  Sin(^diskt)  .  1  [*  t (  \  u\  M  j 

w(t)  — - — --) - /  f(r)  h[w(t  -  t)\  dr 

OJ disk  0) disk  'W'disk  JO 

OJdisk  ~  OJndisk  \J  1 

h  ( t )  =  e~^ndiskt  sin  {ojdisk  t ) 


(3.4.10) 


(3.4.11) 

(3.4.12) 

(3.4.13) 


where 


OJ-n  clamp  — 


^c  lamp 
\  m  clamp 


fcdi 


l^ndisk  — 


disk 


c  = 


Tndisk 
dampdisk 
2  nidisk 


(3.4.14) 

(3.4.15) 

(3.4.16) 


and  v ciamp,  Vdisk  are  the  initial  velocities  of  the  clamp  and  the  disk.  Equations 
(3.4.9)  through  (3.4.11)  give, 


m 


i 


clamp  W clamp  •'O 


[  f(r)sin[w(t  —  T)]dT 
Jo 


OJdisk  m disk 


J ^  f(T)h[w(t-T)]dT^ 


3/2 


(3.4.17) 
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where, 


_  ^ciamp  Sin  {t^nclampi')  ^  C\V  sin  {cOnclampt] 

M n  clamp  ^clamp^n  clamp 


Vdisk  sin(u>diskt) 


^ disk 


(3.4.18) 


The  disk  and  clamp  system  is  similar  to  the  lumped-parameter  structure  we 
discussed  in  section  2.2.1.  Without  much  difficulty,  we  can  show  that  the  Green’s 
function  of  this  system  satisfies  the  conditions  of  theorem  2.4.2  and  2.4.3.  Hence, 
the  Hertz  equation  has  a  unique  solution  and  the  impact  model  is  valid.  The 
numerical  solutions  of  the  impact  model  for  the  actuator  will  be  given  in  the 
next  chapter.  For  simplicity,  this  dissertation  has  dealt  solely  with  the  case  of 
beam.  Extensions  of  the  method  to  multiple  space  dimensions  e.g.  plates,  are 
feasible. 
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Chapter  4 


NUMERICAL  AND  APPROXIMATION  METHODS 


In  chapter  three  we  used  the  Hertz  law  of  impact  to  derive  an  impact  model  when 
impact  involves  a  flexible  structure,  and  we  also  established  the  existence  and 
uniqueness  of  solutions  of  the  Hertz  equation.  Unfortunately,  this  impact  model 
is  nonlinear  and  does  not  admit  a  closed  form  solution.  Hence,  computational 
aspects  must  be  considered. 

Inspired  by  the  contraction  mapping  theorem,  we  develop  a  numerical  method 

using  the  successive  Picard  iterations:  /0,  P/o,  PPfo, . ,  where  the  initial 

condition  /0  is  obtained  from  the  energy  method,  and  P  is  a  contraction  operator. 
Our  experience  is  that  this  method  is  faster  compared  to  the  small-increment 
method,  more  accurate  than  the  energy  method,  and  its  convergence  is  very  fast. 

Although  the  methods  (  small-increment  method  and  energy  method  etc.) 
could  provide  the  numerical  solutions  to  the  Hertz  equation,  these  methods  share 
some  drawbacks,  the  main  one  is  that  parameter  variations  e.g.  varying  initial 
velocities,  demand  a  whole  new  computation  of  the  numerical  solution.  Also,  a 
fairly  large  computational  burden  has  to  be  incurred  for  each  numerical  solution. 

We  propose  a  method  to  alleviate  this  problem.  Let  us  first  observe  from  our 
experience  that  the  impact  period  is  very  short  in  general.  Therefore,  one  may 
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approximate  the  transfer  functions  of  the  systems  to  which  the  impacting  bodies 
belong  by  Taylor  polynomials  of  low  order.  We  explicitly  carry  out  this  compu¬ 
tation  in  the  cases  of  first  and  second  order  Taylor  polynomial  approximations. 
We  show  that  in  the  case  of  the  first  order  approximation,  there  is  a  universal 
ordinary  differential  equation  that  describes  the  impact  behavior  completely  in 
the  sense  that  parameter  variations  only  require  proper  rescaling  of  functions. 
Therefore,  one  can  numerically  solve  this  equation  beforehand,  save  the  results, 
and  can  use  it  to  predict  the  impact  behavior  with  only  a  minimal  computational 
burden.  In  the  case  of  the  second  order  approximation,  there  is  a  two  parameters 
family  of  ordinary  differential  equations  that  govern  the  impact  behavior. 

4.1  A  Numerical  Method 

The  energy  method  is  simple  and  fast,  although  less  accurate  than  the  small- 
increment  method.  Hence,  the  results  obtained  by  the  energy  method  can  be 
used  as  a  good  initial  approximation. 

Inspired  by  the  contraction  mapping  theorem,  we  develop  a  numerical  method 
using  the  successive  Picard  iterations,  /o,  P/o,  PP/o, . ,  where  the  initial  con¬ 

dition  /o  is  obtained  from  the  energy  method,  and  P  is  the  contraction  operator 
defined  by, 

1  rt  3/2 

Pf(t )  =  W0t - -  f(r)L(t  -  r)dr}  ,  Vt  G  [0,T0]. 

TYl  J  0 

Recall  that  the  function  L(t)  is  given  by 

L(t)  =  t  +  mG(x*-,t). 
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The  successive  approximations  to  a  solution  of  (3.2.10)  are  given  by 


1  rt  3/2 

fn(t)  =  K t  -  —  yo  /„_i(r)L(t  -  r)dr]  ,  n-  1, 2,  —  (4.1.1) 

We  take  the  initial  condition  /0  to  be  the  one  given  by  (2.2.18).  Closed  form 
solution  of  the  first  iteration  is, 


fi(t) 


K t  ~  J0  Mr)Ht  ~  T)dT]3/2, 

r  /  4.  Ko  /  ^  1  \ ,  \ 

[”»t  “  n7  A  -  XSmXt) 

_TS  ^s'i'nwkt  ~  WkSinXt. [3/2 

°£ri  }  W*(A2  -  tug)  J 


Vt  €  [0,  Tq  . 


(4.1.2) 


It  turns  out  that  even  this  first  approximation  yields  reasonable  results  in  some 
cases  (see  e.g.  Figure  4.4). 


4.2  Numerical  Examples 

We  consider  some  impact  examples  involving  a  flexible  structure;  a  simply  sup¬ 
ported  beam  on  which  a  spherical  body  impacts  at  the  center.  Solution  of  the 
Hertz  equation  depends  on  several  parameters.  The  following  two  examples  are 
used  to  show  how  the  impact  force  is  affected  by  the  approach  velocity  vq  and 
material  properties,  and  these  results  give  us  some  “feel”  for  the  impact  problem. 
The  model  parameters  used  in  the  numerical  simulations  are  given  in  Table  3.1. 


Remark  4.2.1  The  purpose  of  the  first  example  is  to  address  the  question  of  how 
the  impact  force  is  affected  by  the  approach  velocity.  We  use  different  approach 
velocities  while  keeping  the  material  properties  fixed.  Force  profiles  are  plotted  in 
Figure  f.l. 
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m\.  mass  of  the  impactor  (steel  ball) 

0.2  kg 

v0:  velocity  of  impactor 

0.1  m/s 

l  :  length  of  the  beam 

1.0  m 

h:  thickness  of  the  beam 

0.02  m 

b:  width  of  the  beam 

0.06  m 

Table  4.1:  Model  parameters 


Figure  4.1:  Impact  forces  by  different  velocities;  material:  aluminum 

For  a  given  pair  of  impacting  objects  the  impact  forces  and  impact  durations 
are  determined  by  the  approach  velocity  Vq.  Smaller  v'0s  result  in  smaller  impact 
forces  and  longer  impact  durations,  while  the  larger  v'0s  result  in  greater  impact 
forces  and  shorter  impact  durations.  The  approach  velocity  Vo  is  usually  control¬ 
lable  in  robotics  applications.  Hence,  it  can  be  effectively  exploited  in  the  control 
design  [59]. 

Remark  4.2.2  The  second  example  illustrates  how  the  impact  force  is  affected 
by  material  properties.  We  use  three  different  materials  typically  used  in  robotics. 
The  approach  velocity  is  fixed  and  the  impact  force  profiles  are  plotted.  The  fol- 
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300, 


Figure  4.2:  Impact  forces  using  different  materials;  velocity:  v0  =  0.1  m/s 

lowing  conclusions  can  be  drawn  from  Figure  f.2.  For  a  given  approach  velocity 
Vo,  the  force  magnitudes  and  impact  durations  are  determined  by  material  prop¬ 
erties  of  both  objects.  The  softer  material  (aluminum)  results  in  smaller  impact 
force  and  longer  impact  duration  while  the  harder  one  (steel)  results  in  greater 
impact  force  and  shorter  impact  duration.  This  relationship  leads  to  the  con¬ 
cept  of  passive  compliance.  When  the  approach  velocity  cannot  be  controlled  to 
be  arbitrarily  small,  the  passive  compliance  can  be  introduced  by  the  use  of  soft 
contact  surfaces  (such  as  soft  force  sensors),  thereby  reducing  the  impact  force 
greatly  [47]. 

4.2.1  A  simply  supported  beam 

We  compare  the  three  numerical  methods  discussed  before.  We  consider  a  sim¬ 
ply  supported  beam  case  for  which  parameters  are  given  in  Table  3.1.  Impact 
occurs  at  the  center  of  the  beam.  Figure  4.3  shows  the  solutions  obtained  from 
the  energy  method,  the  small-increment  method  and  three  iterations  using  the 
solution  obtained  from  the  energy  method  as  the  initial  data.  We  note  that  the 
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Figure  4.3:  Central  impact  on  a  simply  supported  beam 

energy  method  yields  a  large  error.  In  order  to  apply  the  Picard  iteration  method 
to  this  case  we  need  first  to  show  that  the  Green’s  function  associated  with  this 
case  is  uniformly  bounded.  Details  are  given  in  Appendix  A.l.  The  result  ob¬ 
tained  from  three  iterations  is  very  good,  and  the  computational  complexity  is 
just  1/3  that  of  the  small-increment  method. 

4.2.2  A  cantilevered  beam 

In  this  example,  we  consider  a  cantilevered  beam  and  the  system  parameters  are 
the  same  as  in  Table  3.1.  Impact  occurs  at  the  tip  of  the  beam.  We  see  that 
fairly  large  errors  occurred  by  using  the  energy  method.  On  the  other  hand, 
the  Piccard  iteration  method  gives  a  good  closed  form  result  even  after  just  one 
iteration.  Figure  4.4  also  shows  the  fast  convergence  of  this  algorithm.  Again, 
the  Green’s  function  for  this  case  is  uniformly  bounded  as  proven  in  Appendix 
A. 2. 
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Figure  4.4:  Tip  impact  on  a  cantilevered  beam 


IIHBH 

Material 

Mass  kg 

Radius  m 

Velocity  m/s 

Body  1 

Phosphor  Bronze 

0.0114 

00 

0.03 

Body  2 

Steel 

1.044 

0.05 

-0.01 

Table  4.2:  Impact  parameters  of  motor 


4.2.3  A  smart  actuator 

We  now  revist  the  example  of  a  smart  actuator  discussed  in  chapter  2.  A  typical 
impact  force  profile  is  given  in  Figure  4.5,  which  is  useful  in  the  design  proce¬ 
dures  such  as  selecting  materials  etc.  Venkataraman  [58]  has  done  extensively 
numerical  simulations,  and  found  the  results  by  small-increment  method  and 
2-steps  Piccard  iteration  method  are  essentially  same  but  the  computation  of 
the  latter  is  just  1/4  of  the  former  method. 

4.3  Approximation  Methods 

Numerical  data  obtained  from  the  impact  model  provides  useful  information 
as  illustrated  above.  These  numerical  methods  share  some  drawbacks,  one  of 
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Figure  4.5:  Impact  force  profile 

which  is  that  each  situation,  e.g.,  resulting  from  varying  initial  velocities,  has  to 
be  numerically  solved  separately.  Also,  a  fairly  large  computational  burden  is 
incurred  for  each  numerical  solution. 

Our  objective  in  this  section  is  to  show  that  it  is  possible  to  take  advantage 
of  the  fact  that  the  impact  period  is  very  short  in  general.  Therefore,  one  may 
approximate  the  transfer  functions  of  the  systems  to  which  the  impacting  bodies 
belong  by  Taylor  polynomials  of  low  order.  We  carry  out  this  computation  ex¬ 
plicitly  in  the  cases  of  first  and  second  order  Taylor  polynomial  approximations. 
We  show  that  in  the  case  of  the  first  order  approximation  there  is  a  univer¬ 
sal  ordinary  differential  equation  that  describes  the  impact  behavior  completely. 
Therefore,  one  can  solve  this  equation  numerically  beforehand,  save  the  results, 
and  use  them  to  predict  the  impact  behavior  with  only  a  minimal  computational 
burden.  In  the  case  of  the  second  order  approximation,  there  is  a  two-parameter 
family  of  ordinary  differential  equations  that  govern  the  impact  behavior. 


57 


In  order  to  motivate  the  problem,  for  the  sake  of  simplicity,  let  us  consider 
two  flexible  bodies  in  linear  motion  as  shown  in  Figure  4.6.  This  could  be  a 
simplified  model  for  a  robot  and  a  work- piece  [43,  60],  where  is  the  mass  of 
the  ith  body,  c;  and  k{  are  the  associated  viscous  frictional  coefficient  and  stiffness 
respectively,  and  u(t)  is  the  control  input  (force).  Let  f(t)  denote  the  impact 
force  between  the  bodies.  The  dynamical  equations  for  the  system  during  the 


Figure  4.6:  A  model  of  robot  and  workpiece 


impact  are: 

m,ix\(t)  +  CiXi(t)  +  kiXi(t)  =  u(t )  -  f(t), 

m2x2(t)  +  c2x2(t)  +  k2x2(t)  =  f(t).  (4.3.1) 

Without  loss  of  generality,  we  assume  that  the  approach  velocity  of  the  first 
body  just  before  contact  is  v0  >  0  and  that  the  second  body  is  initially  at  rest. 
The  Hertz  law  of  impact  is  assumed  valid,  i.e. 

a(t)  =  K\l(t)f'3.  (4.3.2) 

Let  us  assume  that  the  external  forces  acting  on  the  bodies  are  negligible  in 
comparison  to  the  impact  force.  We  can  write  the  displacements  of  the  bodies 
as 

Xi(t)  =  h0(t)  -  f  h^t  -  r)/(r)dr, 

Jo 
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x2(t)  =  [  h2{t  -  r)/(r)rfr. 

jo 


(4.3.3) 


where  h0(t)  is  a  function  of  initial  condition  and  h0(t)  =  v0mihi(t).  The  Green’s 
functions  hi  are  given  by 


hi(t) 


LOi 


k  ym 

where  i  =  1,2  and  and  (j  are  defined  by 


e  ,\f\  -  Q  t),  (4.3.4) 


•4  =~,  -  =  20*- 

rrij  m, 


Hence,  the  relative  approach  can  be  written  as 


a(t) 


ho(t)  -  f  hi(t-T)f(r)dT  -  f  h2(t  —  r)/(r)dr, 

Jo  Jo 

ho(t)  -  [  h{t-  r)/(r)dr,  (4.3.5) 

Jo 


where  h{t)  :=  h\{t)  +  h2{t). 

Equations  (4.3.2)  -  (4.3.5)  lead  to, 

K[f{t)}2/3  =  h0(t)  -  [  h{t  —  r)/(r)dr  (4.3.6) 


If  the  Green’s  function  h{t)  belongs  to  C°°[0,  T]  (this  condition  is  satisfied  in 
many  applications),  then  h(t)  can  be  expanded  in  the  form  of  a  Taylor  series: 


h{t)  —  kit  +  —  k2t2  +  + 


(4.3.7) 


where  the  coefficients  ki  are  determined  by  h{t)  and  are  finite.  Since  the  impact 
duration  is  short,  one  can  frequently  approximate  h{t)  by  its  first  or  second  order 
Taylor  polynomials  very  well. 
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4.3.1  First-order  approximation 


The  solution  of  nonlinear  equation  (4.3.6)  can  be  greatly  simplified  if  we  use 
the  first-order  approximation.  Somewhat  surprisingly,  equation (4. 3. 6)  renders 
a  universal  ordinary  differential  equation  (UODE)  in  this  case.  The  first-order 
approximations  of  functions  ho{t)  and  h(t)  are  v0t  and  k\t  respectively.  By 
plugging  these  approximations  into  equation  (4.3.6)  we  obtain 

Klf(t)}2/3  =  v0t-  f  ki(t-T)f(r)d,T.  (4.3.8) 

Jo 

Because  the  relative  approach  a(t)  =  Kf2^(t),  replacing  f(t)  by  a(t)  in  equation 
(4.3.8)  gives 


a(t)  =  v0t-  [  -  T)a3/2(r)dT. 

Jo  K 

Differentiating  equation  (4.3.9)  twice  leads  to, 


d2a(t ) 
dt 2 


-*l(I)3/2Q3/2(t). 


(4.3.9) 


(4.3.10) 


By  our  assumptions,  the  initial  conditions  are  cx(0)  =  0,  and  d(0)  =  v0.  Equation 
(4.3.10)  can  be  integrated  once  to  obtain 

(^)2-»i i  = -|*i(^)3/2“5/2W-  (4.3.11) 


From  equation  (4.3.11),  the  maximum  relative  approach,  amax,  can  be  imme¬ 
diately  obtained  by  letting  a(t)  =  0,  and  the  maximum  impact  force  is  easily 
obtained  from  (4.3.2): 


fv 


,5*3%2/5 
4  fci  ’ 


[- 
L  A 


4  kxK 


]3/5. 


(4.3.12) 
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Next,  we  show  that  equation  (4.3.11)  can  be  put  into  a  universal  form  by  scaling 
both  the  time  and  magnitude  of  a(t).  Let  us  introduce  two  new  variables  (5  and 
t  by 


/ 3(t )  =  A  a(t),  t  —  nr,  (4.3.13) 

where  (3  and  r  will  play  the  roles  of  the  dependent  and  the  independent  variables, 
A  and  n  are  two  constants  to  be  determined  shortly.  It  is  easy  to  show  that 
=  A Replacing  a  and  t  by  the  expressions  in  (4.3.13)  we  obtain  from 
(4.3.11) 

(^)2-AV»2  =  -iA(4)3/2-L/i(r)5/’.  (4.3.14) 

Let  us  now  choose  A  and  /i  such  that 

AV«o  =  1,  ki(-j^)3/2  /  V\  =  1.  (4.3.15) 

Equation  (4.3.14)  now  takes  the  universal  form 

.  ^  =  v/l  -  j8(r)«/2,  (4.3.16) 

subject  to  the  initial  condition  /?( 0)  =  0.  This  completes  the  construction.  □ 
The  universal  ordinary  differential  equation  (4.3.16)  has  a  unique  solution  with 
the  given  initial  condition,  and  can  be  solved  numerically.  The  result  is  plotted 
in  Figure  4.7.  Once  the  solution  of  this  equation  is  obtained,  any  impact  prob¬ 
lem  can  be  solved  immediately  by  using  the  following  relation.  From  equation 
(4.3.15),  we  can  solve  for  A  and  /r  respectively: 

a  =  1<5>2^1/5'  (43'17) 
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Let  Tmax  be  the  time  before  the  solution  of  (4.3.16)  crosses  the  time  axis  once 
again.  Then  the  impact  duration  T  for  the  original  problem  is 

T  =  HTmax  =  ]1/5rmaz-  (4.3.18) 

The  impact  force  can  be  obtained  from 

/(*)  =  (^)3/2/3(</M)3/2  =  (4-3.19) 

By  scaling  back  the  time  from  r  to  t,  we  can  recover  the  actual  force  f(t).  Let 


Figure  4.7:  Numerical  solution  of  UODE 

us  apply  this  method  to  the  system  shown  in  Figure  4.6  for  which  the  model  and 
impact  parameters  are  given  in  Table  3.2. 

Case  A:  The  two  impacting  bodies  are  assumed  to  be  made  of  steel.  Figure 
4.8  shows  that  the  approximate  solution  is  practically  indistinguishable  from  the 
exact  solution  (numerically  computed  with  the  small  increment  method). 

Remark  4.3.1  The  first-order  approximation  has  an  interesting  physical  inter¬ 
pretation.  Note  that  k\  —  (1/rai  +  l/m2).  This  means  that  we  treat  the  two 
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m\: 

2.0  kg 

Ci  : 

20.0  N.s/m 

k\’ 

104  N/m 

ra2: 

Q.§kg 

c2: 

15  N.s/m 

&2: 

104  N/m 

v0: 

0.1  m/s 

Table  4.3:  Model  parameters  of  flexible  bodies 


Figure  4.8:  Solutions  by  lst  order  approximation  and  small-incre.  methods:  Case 
A 
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Figure  4.9:  Solutions  by  lst  order  approximation  and  small-incre.  methods:  Case 
B 
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flexible  bodies  as  if  they  were  two  single  rigid  bodies  in  free  space.  When  the 
impact  duration  is  relatively  large  (  low-speed  impact,  softer  materials  contact¬ 
ing  etc.),  the  effects  of  other  factors  such  as  damping  cannot  be  ignored.  The 
first-order  approximation  will  introduce  some  errors. 

Case  B:  For  the  purpose  of  comparison,  we  increased  the  damping  coefficients  q 
by  a  factor  of  10.  Results  are  shown  in  Figure  4.9.  The  solution  obtained  from 
the  first-order  approximation  now  begins  to  deviate  from  the  exact  solution. 
This  motivates  us  to  consider  the  development  of  a  second  order  approximation. 

4.3.2  Second-order  approximation 

The  second-order  approximation  of  the  function  h(t)  is  kpt  +  \k2t2,  where  kx 
0,  k2  ^  0,  and  the  approximation  of  hflt )  is  vfl  +  \kfl2.  We  substitute  these 
approximations  into  the  nonlinear  equation  (4.3.6): 

K[f(t)]2/:i  =  v0t  +  ]-k0t2 -[  [kflt-T) +  ]-k2(t-T)2} f{r)dT.  (4.3.20) 

Z  Jo  z 

Using  the  relation  (4.3.2),  replace  /(•)  by  a(-), 

a(t)  =  v0t  +  \kQt2  -  [  {^)3/2[ki(t  -  t) +  ^-k2(t  -  r)2]a3/2(T)dT. 

Z  Jo  K  Z 

Hence  we  obtain 

^  (4.3.21) 

where,  k[  =  kx(^)3^2,  and  k'2  =  k2(j^)3^2. 

Let  us  introduce  two  new  variables: 

/ 3(t )  =  A  aft),  t  =  piT.  (4.3.22) 
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From  equation  (4.3.21)  we  obtain, 


d3P{r) 


=  +  WMA72] 


A/i  dr 


=  -l*4^d^-k'4^3'2  ^ 

Let  us  first  assume  that  k2  >  0,  and  choose  the  parameters  A  and  p  such  that 


h’^-V  k!  —  —  1 

2*vx~  ’ 

Equation  (4.3.23)  now  becomes  a  universal  differential  equation 


d3/3(r) 
dr 3 


(4.3.24) 


(4.3.25) 


with  the  initial  conditions  given  by 


m  =  0; 


d0(r) 


t=o  -  A/ru0  -  71; 


d  ^  lr=o  =  A/i2fc0  =  72.  (4.3.26) 


Hence,  we  have  a  universal  ordinary  differential  equation  (4.3.25),  and  a  two- 
parameter  family  of  initial  data,  parameterized  by  71  and  72. 

Equation  (4.3.24)  can  be  solved  for  p  and  A  respectively: 

3  k\  3kx  .  .3.0  k[4 

"=2J'=2V  X  =  {2>W 


^  2  k'2  2kff  v2'  kq 

/3\4^15  /  ^1*  7, 

71  =  2  72  =  2  fc°- 
ft  2  ^  A/2 


(4.3.27) 


3 

Remark  4.3.2  HTien  k2  <0,  let  k'2-^j  —  —  1.  By  the  same  procedure,  we  obtain 


the  universal  differential  equation 


d3P(r)  ro7Z\dP(T) 


(4.3.28) 


with  the  same  initial  conditions  (4-3.26),  where  the  parameters  are , 


3  k[  3  ki 


,3.3k'i 


—  ___L  —  \  — 

1  2  k2  2  k2  ’  V  &22’ 

/3.4 /cf  /3.4kf 

71  =  72  =  -b)  pi^o- 

4  /u2  ^  A/2 


(4.3.29) 
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Once  again,  the  universal  differential  equation  (4.3.25)  with  two-parameter  7x 
and  72  can  be  solved  numerically.  A  table  can  be  created  to  store  the  results 
with  different  71  and  72,  and  any  impact  problem  can  be  immediately  solved  by 
using  data  from  the  table.  Obviously,  higher-order  approximations  may  increase 
the  accuracy,  but  only  at  the  expense  of  computational  load.  Let  us  apply  the 
second-order  approximation  to  the  above  examples. 

Case  C:  For  the  purpose  of  comparison,  we  increased  the  damping  coefficients  c* 
by  a  factor  of  50  .  It  is  clear  that  k2  =  —  (ci/mf  +  c2/m\)  <  0;  For  the  given 
model  parameters,  a  table  is  generated  with  two  variables  71  and  72.  It  is  seen 
in  Figure  4.10  that  the  result  obtained  by  using  the  second  order  approximation 
agrees  with  true  results  to  a  high  degree  of  accuracy. 

In  summary,  an  approximation  method  has  been  developed  for  the  analysis  of 


Figure  4.10:  Solutions  by  lst  and  2 n<i  order  approximation  and  small-incre.  meth¬ 
ods 
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dynamics  of  flexible  bodies  undergoing  impact.  The  first-order  approximation 
yields  a  special  function  which  can  be  used  for  analytical  and  computational 
purposes.  This  approximation  seems  to  give  results  which  are  very  close  to  the 
exact  solutions.  A  second  order  approximation  was  developed  as  well,  and  it 
leads  to  a  two  parameter  family  of  ordinary  differential  equations  of  which  the 
solutions  contain  universal  features  of  impact  problems. 


68 


Chapter  5 


NONLINEAR  OPTIMAL  CONTROL  OF  IMPACT 
FORCES 


In  this  chapter,  we  study  some  control  issues  of  impact  dynamics.  We  will  study 
a  control  problem  where  a  linear  system  is  subjected  to  a  series  of  impact  forces. 
The  impact  forces  are  treated  as  disturbances  to  the  system  and  modeled  as 
finite  duration  events  using  the  theory  developed  in  chapters  3  and  4.  Among  the 
motivating  factors  is  the  need  to  study  the  control  problems  related  to  mechanical 
systems  subject  to  impact  forces,  e.g.,  active  control  of  the  suspension  system 
of  a  vehicle  against  irregularities  of  the  road,  stabilization  of  an  antenna  on 
the  space  station  subject  to  impact  from  space  debris,  or  active  damping  of 
vibrations  of  flexible  structures  caused  by  impact  forces  [10,  17].  A  reasonable 
control  objective  in  all  these  problems  is  to  design  a  stabilizing  controller  to 
minimize  the  effects  from  the  disturbances  to  the  controlled  outputs  in  some 
sense.  This  problem  in  turn  can  be  studied  in  the  framework  of  control 
theory. 

An  important  paradigm  in  control  synthesis  is  the  control  problem  in¬ 
troduced  by  Zames  [18].  In  this  formulation  the  disturbances  belong  to  a  ball 
in  a  certain  function  space,  and  a  quadratic  cost  function  is  minimized  for  the 
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worst  disturbance  in  this  set.  Various  theories  of  control  problems  for  linear 
systems  have  been  developed  so  far  by  many  researchers  [19,  20,  21,  22,  23]. 
Solutions  to  these  problems  can  be  obtained  via  frequency  domain  methods,  or 
recently  developed  state-space  methods. 

Recently,  several  researchers  have  attempted  to  extend  the  H '<*,  control  the¬ 
ory  to  the  case  of  nonlinear  systems.  Ball  and  Helton  [25,  26].,  from  a  viewpoint 
of  operator  theory,  discussed  H0 0  control  theory  of  nonlinear  systems  for  the 
first  time.  Basar  and  Isidori  [27,  28]  have  connected  the  control  theory  of 
nonlinear  systems  with  dynamic  game  theory.  In  this  setting,  one  is  naturally 
led  to  a  nonlinear  partial  differential  equation  known  as  the  Isaacs  equation.  A 
straightforward  application  of  the  theory  (see  e.g.,  [29])  to  the  case  of  a  plant 
modeled  by  an  affine  nonlinear  system  shows  that  once  a  solution  of  the  ap¬ 
propriate  Isaacs  equation  is  found,  a  (full-information)  feedback  law  providing 
disturbance  attenuation  (in  the  sense  of  the  L2  gain)  can  be  computed  right  away. 
Van  der  Schaft  [30,  31]  analyzed  the  relation  of  the  L2  gain  between  nonlinear 
systems  and  their  linearization,  and  gave  a  sufficient  condition  for  the  existence 
of  smooth  AToo  state  feedback.  In  addition,  using  the  dissipative  system  theory 
[32,  33],  Van  der  Schaft  has  discussed  the  relation  between  the  L2  gain  and  the 
Hamilton- Jacobi  equation  in  [34]. 

The  disturbance  attenuation  of  impact  forces  is  a  nonlinear  control  problem 
due  to  the  nonlinearity  of  the  impact  model,  which  is  different  from  the  prob¬ 
lems  discussed  in  [36,  34,  35].  It  is  not  affine  in  the  disturbance  input,  and 
the  linearization  around  the  equilibrium  does  not  exist.  Hence  this  nonlinear 
problem  cannot  be  solved  by  the  above  proposed  methods,  but  some  analysis 
can  be  carried  over  to  this  problem.  Our  approach  is  based  on  dynamic  game 
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theory  [27,  29].  In  this  setting,  the  control  problem  naturally  becomes  a  mini¬ 
max  optimization  problem  of  a  cost  function  L(u,  v)  with  two  players,  where  u 
is  the  control  and  v  is  the  disturbance.  Roughly  speaking,  the  player  u  tries  to 
minimize  L(u,  v)  over  U,  while  the  player  v  tries  to  maximize  L(u,  v )  in  V  si¬ 
multaneously,  where  U,  V  are  the  constraint  sets  of  u  and  v  respectively.  [27,  29] 
show  that  the  (sub)  optimal  disturbance  attenuation  problem  has  a  solution  for 
a  given  7  >  0  if  the  minimax  optimization  problem  admits  a  saddle  point.  We 
show  that,  due  to  the  nature  of  impact  dynamics,  the  saddle  point  may  not 
exist  in  this  control  problem.  Motivated  by  the  dynamic  game  approach,  if  the 
information  of  the  disturbance  v  is  assumed  to  be  known  a  priori  (e.g.,  sensors 
can  predict  the  impact  velocity  before  impact),  one  can  seek  an  optimal  strategy 
by  using  this  information  such  that  a  certain  attenuation  level  7  is  achieved.  A 
procedure  is  given  to  compute  the  optimal  strategy  u(v),  and  the  optimal  at¬ 
tenuation  level  7*.  Finally,  by  taking  advantage  of  the  fact  that  the  duration  of 
each  impact  force  is  usually  very  short,  we  derive  a  series  of  approximation  mod¬ 
els  for  the  original  nonlinear  system.  We  show  that  the  higher  order  terms  are 
negligible,  hence  a  special  linearization  of  the  nonlinear  system  can  be  obtained 
by  using  the  leading  term  as  an  approximation.  Thus  the  nonlinear  problem  can 
be  approximated  by  a  linear  one.  The  complete  analysis  and  the  solution  of  this 
linear  control  problem  are  discussed  in  detail  in  the  next  chapter. 
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5.1  Dynamic  equation 


Let  l2  denote  the  Hilbert  space  of  square-summable  complex-valued  sequences 
{xk,  k  >  0}  with  the  norm  defined  by 

00 

IMk  :=  (E  M2)1/2-  (5.1.1) 

k= 0 

Let  L2  denotes  the  Hilbert  space  of  all  complex-valued  Lebesque  measurable 
functions  x{t)  defined  on  [0,  00)  with  the  property  that  |x|2  is  Lebesque  inte¬ 
grate.  The  norm  is  defined  by 

r  00 

M\l2  ■■=  (Jo  \\x(t)\\2dt)1/2.  (5.1.2) 

For  simplicity,  we  consider  a  single-input  and  single-output  (SISO)  linear  system 
subject  to  a  series  of  impact  forces,  occurring  at  constant  time  intervals  kT0,  k  = 

0,1,2, . ,  T0  >  0.  The  analysis  and  theories  developed  here  can  be  easily 

extended  to  MIMO  systems  and  some  nonlinear  systems.  We  further  assume 
that  a  state-space  form  of  the  SISO  system  is  given  by 

N 

x(t)  =  Ax(t )  +  Biu(t)  +  B2  E  f(vk,  t  ~  &T0)  (5.1.3) 

k= 0 

where  x  €  Rn  are  state  variables,  A,Bi,B2  are  constant  matrices  with  appro¬ 
priate  dimensions,  N  >  0  is  a  finite  integer,  andf(vk,t  —  kT0 )  is  the  impact  force 
generated  by  the  kth  impact,  which  is  parameterized  by  the  variable  Vk .  We 

assume  that  |ujt|  <  M,k  =  1,2, . iV  for  some  finite  number  M  >  0.  The 

function  f(vk,t  —  kT0 )  is  a  causal  nonlinear  function  derived  from  the  Hertz  law 
of  impact,  which  is  given  by 

K[f(vk,t  -  kT0)]2/3  =  vk(t-kT0)-[  G(t-r)f(vk,r),  t  >  kT0 

JkTo 

f(vk,t-kT0 )  =  0,  Vt  <  kT0.  (5.1.4) 
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where  G(t)  is  the  Green’s  function  of  the  system  and  is  assumed  to  be  known. 
Since  there  is  no  closed-form  solution  to  the  nonlinear  model  (5.1.4)  in  general, 
and  various  numerical  methods  to  solve  this  nonlinear  equation  are  developed 
in  chapter  3.  It  is  easy  to  show  that  if  v  E  l2 ,  then  impact  forces  J2k=o  f(vk >  t  — 
kT0)  £  L2. 

5.2  Formulation  of  the  control  problem 

In  this  subsection,  we  formulate  a  nonlinear  optimal  control  problem  by  speci¬ 
fying  a  control  objective  and  stating  the  general  assumptions. 

5.2.1  The  Z/2  gain  of  a  nonlinear  system 

Let  us  consider  the  system 

N 

x(t)  =  Ax(t)  +  Biu(t)  +  B2  fivk ,  t  -  kT0) 

k= 0 

z(t)  =  Cx(t )  +  Du(t)  (5.2.1) 

where  z(t)  is  the  controlled  output  vector,  and  C,  D  are  constant  matrices  of 
appropriate  dimension.  If  we  treat  the  impact  forces  as  disturbances  to  the 
system  (5.2.1),  a  nonlinear  control  problem  for  disturbance  attenuation  can 
be  formulated.  A  word  concerning  terminology  is  certainly  in  order  here,  because 
the  f/oo  norm  is  defined  as  a  norm  on  transfer  matrices  and  hence  does  not 
directly  generalize  to  nonlinear  systems.  However,  when  translated  to  the  time- 
domain,  the  Hoo  norm  is  none  other  than  the  L2-induced  operator  norm  from  the 
input-output  viewpoint.  This  latter  norm  is  extended  to  the  nonlinear  setting 
by  specifying  that  the  gain  should  not  exceed  a  prespecified  limit. 
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Definition  5.2.1  Let  7  >  0.  System  (5.2.1)  is  said  to  have  L2  gain  less  than 
or  equal  to  7  if 

r\m\\2<'l2'tvl  (5.2.2) 

,'0  k= 0 

for  all  T  >  0  and  all  J2k=ovl  <  °°>  subject  to  the  initial  condition  rr(0)  =  0. 

Now  let  7  >  0  be  a  fixed  number.  Then  in  the  standard  nonlinear  (sub)  op¬ 
timal  control  problem  (for  disturbance  attenuation  level  7)  one  seeks  a  feedback 
controller  u  such  that  the  closed-loop  system  has  L2  gain  less  than  or  equal  to 
7,  i.e.,  such  that  inequality  (5.2.2)  holds.  In  a  dynamic  game  setting,  the  cost 
function  is  defined  as 

Ly(u,v)=  f  \\z(t) ||2  — 72I>i  (5.2.3) 

J0  k= 0 

The  requirement  that  the  L2  gain  does  not  exceed  to  7  is  equivalent  to  existence 
of  a  saddle  point  {u*,v*)  for  the  dynamic  game 

L7(u*,t))  <  L7(u*,  v*)  <  L7(u,  v*), 
and  L-f(u*,v*)  <  0.  (5.2.4) 


5.3  Main  results 

Before  we  study  the  nonlinear  control  problem,  it  is  necessary  to  make  some 

simplifications.  Since  the  impact  force  f(vk,  t  —  kT0 ),  k  =  0, 1, 2, . ,  N  has 

no  closed- form  expression  in  general,  it  is  impossible  to  analyze  and  solve  the 
control  problem  (5.2.1)-(5.2.2)  explicitly. 
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5.3.1  An  approximation  method 

In  chapter  4,  we  developed  first-order  and  second-order  approximations  for  im¬ 
pact  force  functions.  The  first-order  approximation  is  governed  by  a  universal 
differential  equation  and  it  essentially  captures  all  features  of  impact  dynamics. 
We  use  the  first-order  approximation  to  replace  f(vk,t  —  kT0).  To  be  explicit, 
let  us  take  k  =  0.  Now  f(v0,t )  is  given  by  the  equation 

/(« o,t)  =  (l/K)3/2a3/2(v0,t)  (5.3.1) 


where  K  >  0  is  a  constant  which  depends  on  the  local  geometry  of  the  region  of 
contact  and  the  material  properties  of  the  contacting  bodies.  In  general  K  «  1, 
and 


<x(v0,t)  =  jP(t/r), 


(5.3.2) 


where  A  and  /x  are  given  by 


(5.3.3) 


In  (5.3.3),  fci  is  a  constant  and  /?(t)  is  governed  by  the  universal  differential 
equation 

d0(r) 


dr 


=  -  WF*,  (8(0)  =0. 


From  equations  (5.3.2)-(5.3.3),  we  have 


a(v0  ,t)  =  jP{t/n)  =  £Uo /sp(bvl/st), 


(5.3.4) 


(5.3.5) 


where  the  constant  b  is  given  by 


4  k2  , 

* = [(g>2|yi/5- 


(5.3.6) 
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Since  K  «  1,  it  follows  that  b  »  1.  Here  the  constant  b  can  be  thought  of 
as  an  indicator  of  the  time  scale  of  impact  dynamics.  By  plugging  a(vo ,t)  into 
equation  (5.3.1),  we  finally  arrive  at 

f{v o,  t)  =  (l/K)3/2a3/2(v0,  t)  =  avo/5r(bvl/5t)  (5.3.7) 

where  r(-)  =  /?3/2(-)  is  called  the  normalized  force.  Numerical  solutions  of  (5(r) 
and  r(r)  are  given  in  Figure  5.1.  By  using  the  first-order  approximation,  the 


Figure  5.1:  Numerical  solutions  of  {3(t )  and  r(r) 
nonlinear  system  (5.2.1)  becomes 

N 

x{t)  =  Ax(t)  +  Blu{t)  +  B2'£/avl/5r(bv1k/b(t-kT0)) 

fc=0 

z(t)  =  Cx{t)  +  Du(t).  (5.3.8) 

The  equation  (5.3.8)  has  a  simpler  form  than  (5.2.1),  but  it  is  still  too  compli¬ 
cated  to  analyze.  Next,  we  show  that  it  is  possible  to  take  advantage  of  the 
fact  that  the  impact  period  is  very  short  in  general.  Therefore,  one  may  further 
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approximate  the  nonlinear  system  (5.3.8)  by  Taylor  polynomials  of  low  order.  In 
order  to  motivate  our  approach,  let  us  consider  a  simple  example.  For  a  scalar 
linear  system  subjected  to  a  single  impact  with  velocity  v0,  we  obtain 


x(t)  —  a\x(t)  +  biv^5r(bvl^t),  x(0)  =  0. 


(5.3.9) 


where  ai,bi  are  constants.  x(t)  can  be  solved  explicitly  to  obtain, 

x(t)  =  eait  [  e~aiTbiV^5r(bvl^T)dT,  t  >  0. 

Jo 

Let  us  perform  time-scaling  by  defining 


bvlJbr  —  a. 


Then  x(t)  in  (5.3.10)  becomes 


rbv^5t - iVs0,  K/c  1 

x(t)  =  eait  e  bvo  b^l1  —^r (a) d< 


A.  rbvlJ^t - a^>  g*  <T 

=  eait—Vo  [  e  bvo  r(a)da. 
b  Jo 


a\ 


Expanding  e  bvo  in  the  form  of  a  Taylor  series,  we  have 


— iVs*7 
e  =  1 


a  1  1  a?  o 


(5.3.10) 


(5.3.11) 


(5.3.12) 


(5.3.13) 


k1/5  2!  (bvl/5y 

where  a\  is  the  time  constant  for  the  system  (5.3.9).  A  series  of  approximate 
models  is  given  as 


x(t) 


h.  rbvl,5t 

=  r(ff)[1-^ff  + . ]d° 

=  eaitbiv06i(bvl/5t)  -  eaitb2VQ/b02(bvl/5t)  + . ,  (5.3.14) 


where  0i(t)  and  62  (t)  are  given  by 


Figure  5.2:  Plot  of  nonlinear  function  9\(t) 

Since  impact  dynamics  is  much  faster  than  system  dynamics  in  general,  we  have 
^  <<  1.  Thus  we  can  often  approximate  x(t)  by  low  order  terms  very  well. 
Let  Tmax  be  the  time  before  the  solution  of  equation  (5.3.4)  crosses  the  time  axis 
once  again  (see  Figure  5.1)  and  Tmax  —  Tmax/(bv^5),  define  the  constants  9\  02 
as 

_  fTmax  _  fTmax 

&i=  r{a)da,  92 =  /  crr(a)da  (5.3.16) 

Jo  Jo 

A  plot  of  9\{bvlJbt)  is  given  in  Figure  5.2. 

Remark  5.3.1  Since  0  <  Tmax  «  1,  the  slope  of  9\{hv^bt)  is  very  sharp. 
Notice  that  v09i(bv^5t)  =  0,  Vt  >  0  if  v 0  =  0.  We  can  approximate  the  function 
Vo0i(bv^5t)  by  vq9\  very  well. 
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5.3.2  A  nonlinear  control  problem 

We  now  consider  the  following  infinite  horizon  nonlinear  problem  with  zero  initial 
condition: 

N 

x(t)  =  Ax(t)  +  Biu{t)  +  B2  Y!  avk^r(Pvk/5(t  ~  kTo)) 

k=0 

z(t )  =  Cx(t)  +  Du(t),  (5.3.17) 

where  the  cost  function 

rOO 

J1(u,v)=  \\z\\ldt  -  (5.3.18) 

J  0 

for  a  given  7  >  0.  For  simplicity,  let  w(t)  —  Ytk=oavk^r^)Vkb{i  ~  k^o))- 

In  view  of  the  definition  of  L2  gain,  for  a  given  7  >  0,  we  want  to  find  a 
feedback  control  u  such  that  J7(u,u)  <  0,  Vw  G  V.  It  is  easy  to  show  [27,  29] 
that  this  condition  is  equivalent  to  the  fact  that  a  dynamic  game  admits  the 
following  saddle  point, 

L7(u»  <  L7(u*,u*)  <  L7{u,v*), 
and  L7(u*,v*)  <  0.  (5.3.19) 

In  this  control  problem,  the  saddle  point  (u*,v*)  of  the  dynamic  game  may  not 
exist  due  to  the  nonlinearity  of  impact  dynamics,  hence  it  cannot  be  solved  by 
the  standard  approaches  [27,  29].  Motivated  by  the  condition  (5.3.19),  we  con¬ 
sider  a  related  optimal  control  problem. 

Impact  control  problem: 

Let  us  assume  that  the  disturbance  v  is  known  a  priori.  We  want  to  use  this 
information  to  design  an  optimal  strategy  u(v)  such  that  a  certain  disturbance 
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attenuation  level  7  >  0  can  be  achieved,  i.e.,  the  following  inequality  holds  for 


all  v  eV, 


L7(u(v),v)  <0.  (5.3.20) 

A  sample  problem  of  this  type  is  to  stabilize  an  antenna  on  a  space  station 
subject  to  impact  from  space  debris,  where  we  assume  that  the  debris  coming 
within  a  certain  distance  is  measurable  and  are  modeled  as  a  sequence  of  events. 
Towards  this  end,  we  consider  the  following  max-min  optimization  problem 

max  min  JJu.v)  (5.3.21) 

vev  ueu  v  ' 

under  the  constraint  (5.3.17). 

Remark  5.3.2  The  max-min  optimization  problem  is  much  easier  to  solve  than 
the  minimax  problem.  For  a  given  v,  the  minimization  problem  with  respect 
to  u  can  be  solved  by  the  standard  optimal  control  theory.  The  cost  function 
J7(u(v),v)  can  be  explicitly  evaluated  by  using  the  above  approximation  method. 
The  cost  function  J7(u(v),v)  is  a  nonlinear  function  of  v.  Since  v  €  V  and  V  is 
a  compact  set,  the  maximization  problem  with  respect  to  v  always  has  a  solution. 

We  assume  that  the  state  variables  are  available  and  the  disturbance  v  is  known 
a  priori.  The  following  standard  assumptions  are  made: 

Al)  (A,  Bi)  and  (A,Bf)  are  stabilizable,  and  (C,  A)  is  detectable. 

A2)  D'[C  D]  =  [0  /]. 

Theorem  5.3.1  Let  the  pair  (u,v)  be  a  solution  of  the  max-min  problem 
(5.3.21)  under  the  constraint  (5.3.17).  If  the  cost  function  J^(u,v)  <  0  for  a 
given  7  >  0,  then  the  inequality  (5.3.20)  holds  for  all  v  E  V  and  the  optimal 
strategy  is  given  by  the  solution  of  the  minimization  problem  minu  J7(u,  v). 
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Proof:  Let  u(v)  denote  the  solution  of  the  minimization  problem  of  J7(u,  v),  so 
that  v  becomes  the  only  variable  of  the  function  J^(u{v),  v).  Let  v  be  the  solution 
of  the  maximization  problem  of  J-f(u(v),v)  under  constraint  V .  Then  by  the 
assumption,  we  have  the  following  inequality 


J^(u(v),v)  <  J^(ii,v)  <  0 


(5.3.22) 


where  u  =  u(v). 


□ 


Now  we  discuss  the  problem  of  computing  the  optimal  control  strategy  u(v).  For 
a  given  v,  the  Hamiltanion  of  (5.3.17)-(5.3.18)  is  given  by 


H(x,  u,p )  =  x'C'Cx  +  u'u)  +  p\Ax  +  B\u  +  B2w). 
z 

The  minimum  principle  now  yields 

dH 


du 


=  0,  =$■  u*  =  —B[p 


The  state  and  adjoint  equations  are  given  by 

dH(x,  u*,p ) 


d/dt 


x  — 


P  =  ~ 


X 


dp 

dH{x ,  u*,p) 
dx 


—  Ax  —  B\B[p  +  B2w 
-  -A'p  -  C'Cx 


P 


A  -B1B[ 
-CjCj  -A! 


X 

b2 

+ 

p 

0 

w 


=  H 


x 

P 


+  Pw,  t  >  0 


(5.3.23) 


(5.3.24) 


(5.3.25) 


(5.3.26) 


Note  for  t  >  (N  +  1  )T0,  w(t)  =  0,  and  we  have  the  2 n  x  2 n  Hamiltanion  system 


d/dt 


x 


P 


A  -BXB[ 
-C[C1  —A! 


x 


P 


,  t  >  (N  +  1)T0  (5.3.27) 
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The  algebraic  Riccati  equation  associated  with  (5.3.27)  is 


A'Q  +  QA-Q'B1B[Q  +  C'C  =  0.  (5.3.28) 

The  Riccati  equation  admits  a  unique  solution  Q  =  Q'  >  0  by  assumption  Al) 
and  the  LQG  theory.  It  is  easy  to  show  that  p{{N  +  1)T0)  =  Qx((N  +  1)T0). 
The  linear  system  (5.3.26)  can  be  solved  in  backward  time  with  initial  condition 
[x((N  +  1)T0)  p((N  +  1)T0)]', 

x(t) 

_  P(t) 

+em  T  e~HTPw(T)dT,  0  <  t  <  (N  +  1)T0. 

J(N+1)T0 

We  then  apply  the  approximation  method  developed  in  section  4.3.1  to  simplify 

- \jsHt 

expression  (5.3.29).  A  Taylor  expansion  of  the  function  e  hv*  is  given  as 
— \tkHt  1  1  1 

e  *V  =  / - utHt  +  - — + .  (5.3.30) 

The  following  integral  can  be  written  as  a  series, 
rt  n 

/  e~Hr  Pw(r)dr  =  J3(P6iv*(?i(6wJ/5(t  -  kT0)) 

■*0  k= o 

-HPb2vt/592(bvl/5(t  -  kTQ ))  + . ),  (5.3.31) 


oH(t-{N+l)T0) 


x((N  +  1)T0) 
p((N+l)T0) 


(5.3.29) 


and 


/•(iV+l)To  iT  _  .  .  _ 

/  e~Hr Pw{r)dT  =  £(PW r  -  HPb2vAJ%  + . ),  (5.3.32) 

■'0  l _ n 


N 

I 

fc=0 


Note  that 


/t  rr  /*£  tt  r(iV+l)To 

^  e_  TPw(r)dT  =  J  e~HTPw(T)dr  —  J  e~Hr Pw{r)dT (5.3.33) 

We  can  explicitly  express  x((jV  +  1)T0)  as  a  function  of  v  by  using  equations 
(5.3.29),  (5.3.32)  and  the  initial  condition  x(0)  =  0.  Once  an  expression  for 
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x((N  +  1)T0)  is  obtained,  x(t)  and  p(t )  can  be  expressed  as  functions  of  v  during 
the  interval  [0,  (N  +  1)T0].  For  simplicity,  let 

x{t)  =  fi(t,v),  p(t)  =  f2(t,v),  0  <t<  (N  +  1)T0.  (5.3.34) 


where  /i  and  /2  are  two  nonlinear  functions  of  v  in  general.  Thus  the  optimal 
strategy  u(v)  is  given  by 


u(i) 


-B[p(t) 


f  -B[f2(t,  v),  0<t<(N  +  l)T0 


—B[Qx(t),  t  >  (N  +  1)T0 


(5.3.35) 


By  LQG  theory,  the  control  strategy  u(t)  asymptotically  stabilizes  the  closed- 
loop  system. 

Once  the  expressions  for  p(t),  x(t),  and  x((N  +l)T0)  over  the  interval  [0,  (iV+ 
l)To]  are  obtained  ,  the  cost  function  of  /0°°  can  be  easily  evaluated  as 

roo  r(N+l)To  „  „  roo 

I  \\z\\l dt=  I  \\Cx\\l  +  \\u\\ldt+  I  \\Cx\\l  +  \\u\\ldt  (5.3.3Q) 

Jo  Jo  J(N+l)To 

Let  us  denote  the  first  integral  as  /o7V+1^T°  HCrrUl  +  IHI^di  =  ^i(u)  >  0,  and  from 
the  LQG  theory  the  second  integral  can  be  expressed  as  x((N  +  1  )T0)'Qx((N  + 
1)T0)  =  g2(v)  >  0.  Now  the  total  cost  function  becomes 


J7(u,u(u))  =  gi(v)  +  g2(v)  -  72\\v\\l-  (5.3.37) 


The  maximization  problem  max„  J7(u,u(u))  can  be  solved  for  a  given  7  under 
constraint  V.  Observe  that  the  function  J7(i),w)  decreases  monotonically  as  7 
increases,  hence  we  can  always  find  a  7  such  that  <  0.  The  minimal  7* 

such  that  J7(v,u)  <  0  is  called  the  optimal  disturbance  attenuation  level. 
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In  order  to  illustrate  the  procedure,  let  us  take  a  scalar  linear  system  sub¬ 
jected  to  a  single  impact: 

x(t)  —  2  x(t)  +  u(t)  +  av^5r(bv^5t),  (5.3.38) 

where  v0  is  a  variable  v0  £  [0,  0.1].  The  cost  function  is  defined  as 

r  oo 

J-y(u,  u0)  =  /  (x  +  u)2dt  -  j2Vq.  (5.3.39) 

Jo 

Now  we  follow  the  above  procedures  to  compute  the  optimal  strategy  u(v)  and 
optimal  attenuation  level  7*,  we  have 

H(x,u,p )  =  ]-{x  +  u)2  +  p[2x(t)  +  u(t)  +  avQ/5r(bvl/5t)} 
z 

dH  n _  * 

-7—  =  0  =»  u  —  —  p  —  x 
ou 

H(x,  u*,p)  =  p 2  +  p[x(t)  -  p  +  avl/br(bvl^t)} 


X 


X  =  ^  =  x-p  +  aVo/5r(bvl/5t), 


V  = 


dH 


dx 


=  -p. 


(5.3.40) 


Without  loss  of  any  generality,  let  T0  =  1.  It  is  easy  to  obtain  p(  1)  =  2x(l)  from 
equations  (5.3.40).  Thus 

p(t)  =  2e-(‘-1)a;(l), 

x(t)  —  e1_ta:(l)  +  et~T avybr(bvlQbT),  0  <  t  <  1.  (5.3.41) 


Since  x(0)  =  0,  we  have 


a:(l)  =  e  1  /  e  Tav^5r(bv^5T). 
Jo 


(5.3.42) 


First-order  and  second-order  approximations  are  considered  here  for  the  purpose 
of  comparison.  The  first-order  approximation  yields 

*i(l)  = 


ta 


Pi(t)  =  2e  t-9iv0. 


(5.3.43) 
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a 

b 

9918.555 

7934.8442 

aluminum 

308.8 

247.044 

rubber 

20.7873 

16.63 

Table  5.1:  Impact  parameters 


and  the  total  cost  function  is 

roo  rP  - 

J.}(u,Uo)  =  ^  (x  +  u^dt-^vl  =  2— e\vl  -72Wo-  (5.3.44) 

The  second-order  approximation  yields 

£2(1)  =  e-1(^!?;o  -  j£02v} /5 
p2(t)  =  2e~t(^01vo  -  ^#2Vo/5) 

Jy(u,v 0)  =  2(^1%  -  ^2«o/5)2  -  72wo-  (5.3.45) 

The  parameters  a,  b  are  determined  by  the  local  geometry  of  the  region  of  contact 
and  the  material  properties  of  contacting  bodies.  In  order  to  show  that  the 
approximation  is  valid  for  a  large  range  of  applications,  we  select  three  materials: 
steel,  aluminum  and  rubber.  They  are  representatives  of  a  hard,  medium  and 
soft  material.  The  parameters  are  given  in  Table  5.1. 

Casel:  Steel. 

Jy(u,  v0)  =  9.68v2  -  j2Vq  (5.3.46) 

J2(u,  v0)  =  9.68v2  -  1.6  x  1(TV/5  +  3.2  x  l(nV/5  -  (5.3.47) 

Let  7*  be  defined  as  the  optimal  attenuation  level.  The  equation  (5.3.46)  rep¬ 
resents  the  leading  term  approximation  and  equation  (5.3.47)  represents  the 
second-order  approximation.  The  optimal  attenuation  levels  of  (5.3.46)  and 
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(5.3.47)  are  computed  as  7*  =  3.111,72  =  3.111  respectively. 

Case2:  Aluminum. 

J^(u,v 0)  =  9.68u2  -  72t»o  (5.3.48) 

J2(u,  v0)  =  9.68v2  -  0.117v9/5  +  3.3  x  10~V/5  -  72v2  (5.3.49) 

The  optimal  attenuation  levels  of  (5.3.48)  and  (5.3.49)  are  computed  as  7*  = 
3.111,72  =  3.08  respectively. 

Case3:  Rubber. 


J^(u,v0)  =  9.68v2  — 72Vq  (5.3.50) 

J2(u,v0)  =  9.68v2  -  1.78v9/5  +  0.054v8/5  -  72v2  (5.3.51) 

The  optimal  attenuation  levels  of  (5.3.50)  and  (5.3.51)  are  7J  =  3.111, 7^  =  2.641 
respectively. 

Remark  5.3.3  For  the  hard  and  medium  impact  objects,  the  accuracy  of  the 
leading  term  approximation  is  almost  the  same  as  for  the  second- order  approxi¬ 
mation,  and  the  accuracy  of  the  second-order  approximation  is  improved  by  ap¬ 
proximately  15%  for  soft  impact  materials.  Hence  the  following  conclusion  can  be 
drawn:  For  impact  problems  involving  hard  and  medium  materials  which  consists 
of  a  large  part  of  the  applications  we  discussed  before,  the  leading  term  of  the 
Taylor  series  provides  a  very  good  approximation.  From  the  remark  4-3.1,  this 
approximation  basically  renders  a  linear  system.  For  impact  problems  involving 
soft  materials  and  low  velocity,  we  need  to  keep  the  higher  order  term  in  order 
to  obtain  better  approximations  and  need  to  solve  the  nonlinear  problem  directly. 

In  the  cases  where  the  leading  term  is  sufficient,  the  impact  control  problem  is 
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Definition  6.3.2  [62]  An  eigenvalue  A  o/^(T,  0)  is  said  to  be  ( A(t),B(t ))  un¬ 
controllable  if  there  exists  an  n-dimensional  vector  r)  ^  0  such  that 

^(T,  0)'r)  =  \r]  (6.3.5) 

B'q-'it, 0)'r)  =  0,  Vt€[0,T].  (6.3.6) 


Remark  6.3.3  The  pair  (A(t),  B(t))  is  said  to  be  stabilizable  if  and  only  if  all 
the  eigenvalues  A  o/^(T,  0)  such  that  |A|  >  1  are  (A(t),  B(t))  controllable. 


Let  $(i,  t)  defined  in  (6.2.11)  be  partitioned  into  four  n  x  n  matrices  as  follows, 


$(f,r)  = 


$n(*,F) 

$21  {t,r)  $22  (t,r) 
Obviously  $(t,  r)  has  the  property  for  all  t  and  r, 


(6.3.7) 


$(t  +  T,r  +  T)  =  $(t,r). 


(6.3.8) 


Lemma  6.3.1  The  eigenvalues  of$(t  +  T,  t)  are  independent  oft. 
Proof:  Without  loss  generality,  we  assume  t  E  (0,T).  By  definition, 


$(t  +  T,  t ) 


$H(t  +  T,T+)F$H(T,t) 
<t>H(t,0+)F$H(T,t) 
$-H\T,t)$H(T,0+)F$H(T,t) 
^1(T,t)$(T,0)$ff(T,t). 


Thus,  the  eigenvalues  of  $(t  +  T,t)  are  independent  of  t.  □ 

The  hybrid  state  transition  matrix  $(t  +  T,  t )  displays  same  properties  familiar 
with  periodic  systems.  Hence  we  may  expect  that  there  exists  a  periodic  solu¬ 
tion  of  equations  (6.2.1)-(6.2.2).  For  simplicity,  we  will  assume  that  the  matrix 
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$(T,  0)  has  distinct  eigenvalues.  Let  A  be  the  corresponding  diagonal  Jordan 
form  such  that 


S-^O^T,  0)5(0)  =  A, 


(6.3.9) 


where  the  matrix  5(0)  is  composed  of  eigenvectors  of  <3>(T,  0).  It  should  be  noted 
that  the  derivation  of  results  can  be  applied  to  the  case  of  multiple  eigenvalues 
with  minor  modifications.  Now  partition  5(0)  into  four  n  x  n  matrices 


S(0)  = 


y(o)  no) 

x(o>  n») 

and  partition  A  similarly.  Then  equation  (6.3.9)  can  be  written  as, 


(6.3.10) 


®u(T,0)  $i2(T,  0) 

Y(0)  V(0) 

F(0)  V(0) 

Ai  0 

^21(^0)  ®22{T,  0) 

X(0)  U( 0) 

_X(0)  U(O) 

0  a2 

(6.3.11) 


where  Ai  is  an  n  x  n  diagonal  matrix, 

Ai  =  diag{ \u . ,An} 


(6.3.12) 


and  { Ai , . ,  A„}  are  n  of  the  2n  eigenvalues  of  $(T,  0).  Let  us  define  a  2n  x  2n 

matrix  S(t)  for  t  >  0  by 

5(t)  =  $(t,  0)5(0).  (6.3.13) 

Again,  partition  S(t)  as  follows, 

'  Y(t)  V(t) 

X (<)  U(t) 

Lemma  6.3.2  The  equation  $(t  +  T,t)S(t )  =  S(t) A  holds  for  all  t  >  0. 


S(t)  = 


(6.3.14) 
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Proof; 


Thus, 


$(t  +  T,t)S{t)  =  $(f  +  T,T)$(T,t)5(i) 

=  $(i,0)$(T,f)$(f,0)5(0) 
=  $(i,0)$(T,  0)5(0) 

=  $(t,0)5(0)A 
=  S(t)  A 

<f>(t  +  T,  *)5(f)  =  5(t)A,  Vt  >  0. 


(6.3.15) 

□ 


It  is  easy  to  show  that  the  following  equation  also  holds, 
dS(t) 


dt 


=  HS(t),  t  ±  iT , 


S(t+)  =  FS(t),  t  =  iT,i  =  0,1,2, 


Theorem  6.3.1  Zet  A  6e  an  eigenvalue  of$(t,  t)  and 


y 

x 


(6.3.16) 


be  the  correspond¬ 


ing  eigenvector.  Then  A  1  is  also  an  eigenvalue  o/ <3>(f,  r)  and 

eigenvector  of$(t,  r)'  associated  with  A-1. 

Proof;  Define  a  2n  x  2n  symplectic  matrix 

0  / 


—a; 

v 


is  the 


J  = 


-I  0 


(6.3.17) 


Since  H  is  a  Hamiltonian  matrix  and  F  is  a  symplectic  matrix,  we  have  the 
following  equations, 


H'J  +  JH  0, 
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F'JF 


J. 


(6.3.18) 


On  the  other  hand, 


d 


=  $(t,T)'[H'J  +  JH}$(t,T)  =  0. 


(6.3.19) 


The  possible  initial  conditions  are 


$(r,  r)'J$(r,  t)  =  J,  t  ^  iT, 
<S>(t+,t)'J$(t+,t)  -  F'JF  =  J,  r  =  iT 


(6.3.20) 


Hence, 


<$(£,  r)' J$(t,  r)  =  J,  Vt,r  (6.3.21) 

Equation  (6.3.21)shows  that  $(t,r)  is  non-singular  for  all  t  and  r  implying  that 
none  of  the  eigenvalues  of$(t,r)  vanishes.  Let, 


$(i,r) 

y 

y 

=  A 

X 

X 

then,  from  equation  (6.3.21),  we  have 


*- Ht,r y 

—x 

=  J*(t,T ) 

y 

=  A  J 

y 

=  A 

—  X 

y 

X 

X 

y 

Hence, 


(6.3.22) 


(6.3.23) 


—x 

—x 

*(*,r)' 

y 

=  A"1 

y 

Therefore,  if  X  is  an  eigenvalue  of^(t,r),  then  so  is  A  1.  □ 
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6.3.2  Parameterization  of  periodic  solutions  with  jump 

The  following  theorem  gives  a  parameterization  of  all  periodic  solutions  of  equa¬ 
tions  (6.2.1)-(6.2.2)  in  terms  of  X(t)  and  Y(t)  defined  in  equation  (6.3.14). 

Theorem  6.3.2  Suppose  Y(t)  defined  in  equation  (6.3.14)  non- singular  for 

all  t  G  [0,T].  Then  P(t )  given  by 

P(t)  =  X(t)Y{t)~\  t>  0  (6.3.24) 

is  aperiodic  solution  of  equation  (6.2.1)  and  (6.2.2). 

Proof. 

S(t  +  T)  =  $(i  +  T,  0)5(0) 

=  $(t  +  T,t)$(t,0)S(0) 

=  S(t)  A. 

Therefore, 

Y(t  +  T ) 

X(t  +  T) 

P(t)  =  XitjYit)-1  -  XMYiQ-'YMYit)-1 

=  [-C&Yit)  -  A,X(t)]Y(tyl  -  X(t)Y(tyl 
[AY(t)  -  BiB[X(t)]Y{t)~l 
=  -A! P{t)  -  P(t)A  +  P(t)BiB[P(t)  -  C[C\ 

P{  0+)  =  X(0+)F(0+)"1 

=  x(o)[y(o)-'Y-2b2b'2x(o)]-1 

=  P(0)[I --f-2B2B'2P{0)}-1  (6.3.26) 


Y(t) 

X{t) 


Ai,  t  >  0. 


(6.3.25) 
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It  is  easy  to  show  that  P( 0+)  satisfies  equation  (6.2.2).  Thus  P(t)  satisfies 
equations  (6.2.1)  and  (6.2.2).  It  is  periodic  because  from  (6.3.25), 

P(t  +  T)  =  X{t  +  T)Y(t  +  T)-1  =  =  P(t).  (6.3.27) 

□ 


6.3.3  Analysis  of  periodic  solutions  with  jumps 

The  following  Lemmas(4.3)-(4.6)  characterize  the  properties  of  a  solution  P(t) 


that  depend  on  the  choice  of  n  eigenvalues  Ai, . An  for  Ai.  Let 

Q(i)  =  Y(t)*X(t),  (6.3.28) 

fl(t)  =  Y(t)'X(t),  (6.3.29) 


where  *  is  used  to  denote  the  complex  conjugate  transpose. 


Lemma  6.3.3  If  X*  ^  A j1  for  all  i,j ,  1  <  i,j  <  n,  then  Q(t)  is  Hermitian  and 
Q(t)  is  symmetric. 

Proof:  Let  Xi(t)  and  y fit)  be  the  ith  columns  of  X(t)  andY(t)  respectively,  and 
also  let 


zfit) 


Vi(t) 

Xi(t ) 


The  ij  —  element  of  Q(t)  is  then  described  as 


“ij(t)  =  VittM*) 


(6.3.30) 


and  hence 


UJij  (t)  -  =  Z*(t)J Zj  (t) 


(6.3.31) 
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Noting  that  A*  —  Xj 1  ^  0,  equations  (6.3.15)  and  (6.3.21)  leads  to, 

cjijit)  -  uj'iit)  =  (A*  -  (t)Jzj(t)  -  Xfz*  (; t)Jzj(t )] 

=  (A*  -  A;1)-1*? (*)[*(*  +  T,t)'J-  J^\t  +  T:t)}Zj(t) 

=  0.  (6.3.32) 

The  proof  that  Cl(t)  is  symmetric  is  analogous  to  above.  □ 

Lemma  6.3.4  Assume  that  |Y"(£)|  ^  0  for  all  t  6  [0, T\.  If  A*  ^  A j1  and 
A  i  7^  A  j1  for  all  i,j ,  1  <  i,j  <  n,  then  P(t)  defined  by  equation  (6.3.24)  rea l- 
Proof:  Under  the  assumptions,  P(t)  can  be  rewritten  as, 

p(t)  =  x{t)Y{t)~l 

=  (Yitr'yawit)-1 

=  (Y(t)~1yCl(t)Y  (t)~l .  (6.3.33) 

By  Lemma  (6.3.3),  Ut(t)  and  Q,(t)  are  Hermitian  and  symmetric  respectively. 
Equation  (6.3.33)  then  shows  that  P(t)  is  Hermitian  as  well  as  symmetric.  Thus 

P(t)  is  a  real  matrix.  □ 


Lemma  6.3.5  Suppose  that  (A,Bi)  is  controllable  and  ( C\,A )  is  observable.  If 

| Aj|  <  1  for  all  i  =  1,2, . n,  then  Ll(t)  is  positive  definite. 

Proof:  Let  k  =  0, 1, 2, . 


U(k) 


then  from  equation  (6.3.15) 


Y(t) 


A* 


X{t) 


(6.3.34) 


U(k  +  1)  =  $(t  +  T,t)U(k). 


(6.3.35) 


Let  us  define  a  2n  x  2n  matrix  V  as, 


V  = 


0  0 

-I  0 


Q(t)  can  then  be  expressed  as 


fl(t)  =  -U*(0)VU(0). 


Since  |A*|  <  1  for  all  i  —  1,2, . n,  then  U(k)  — >  0  as  k  ^  oo. 

OO 

Q(t)  =  Yu*(k  +  1)vu(k  +  1)~u*(k)vu(k) 

k= 0 

OO 

=  Y,  u*(k)[${t  +  T,  t)'V$(t  +  T,t)  -  V]U{k) 

k= 0 

Define  a  matrix  M(t,  r)  by, 


M(t,r)  =  $(t,7-)'V$(t,r)  -  V. 


Then 


where 


Since 


=  $(t,T)'[H'V  +  VH}<f>(t,T ) 
=  $(t,T)'N'N$(t,T) 


N 


Ci  0 
0  B[ 


M(r,  r)  =  0,  r  iT, 


M(r+ ,  r) 


F'VF  -  V  = 


0  0 
o  7-2b2b^ 


>  0,  r  = 


(6.3.36) 

(6.3.37) 
Therefore, 

(6.3.38) 

(6.3.39) 

(6.3.40) 

(6.3.41) 

(6.3.42) 
iT.  (6.3.43) 
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Thus, 


M(t,T )  =  J*$(s,T)'N'N<f>(s,T)ds,  r^iT,  (6.3.44) 

M(t,r )  =  J  $(s,  t)'N'N^{s,  r)ds  +  M(t+,  t),  r  —  iT.  (6.3.45) 

Since  ( A,B i)  is  controllable  and  ( C\,A )  is  observable,  it  is  easy  to  show  that 
M(t,r )  >  0  for  all  t  >  t.  Thus  from  equations  (6.3.38)  and  (6.3.39),  it  can  be 
concluded  that  Ll(t)  >0.  □ 

Theorem  6.3.3  Assume  that  for  all  eigenvalues  o/ 4>(T,  0),  |A,|  ^  1,  i  — 

1,2 . 2 n  and  that  |T(£)|  7^  0  for  all  t  E  [0,T].  If  n  eigenvalues  in  Ai  are 

chosen  such  that  |A*|  <  1,  i  =  1, 2 . n,  then  a  periodic  solution  P(t)  given  by 

equation  (6.3.24)  rea(  symmetric  and  positive  definite. 

Proof:  The  assumption  |A*|  <  1  for  alii  =  1,2 . n  guarantees  that  A*  7^  A  j1 

and  Xi  7^  A  J1  for  all  i,j ,  1  <  i,j  <  n.  Therefore  P(t)  is  real  and  symmetric  by 
Lemmas  (6.3.3)  and  (6.3.4)  respectively.  That  P(t)  is  positive  definite  follows 
from  Lemma  (6.3.5)  and  equation  (6.3.34).  1=1 

6.3.4  The  periodic  Lyapunov  equation 

Our  approach  for  the  analysis  of  the  closed-loop  system  (6.2.12)  is  based  on  the 
Lyapunov  method.  As  is  well-known,  the  Lyapunov  equation  plays  an  important 
role  in  the  analysis  of  Riccati  equations  and  the  associated  closed-loop  systems 
[63,  69,  70].  Here  we  extend  some  useful  results  on  the  periodic  Lyapunov  equa¬ 
tion  to  our  particular  problem,  namely,  to  the  periodic  Lyapunov  equation  with 
jumps. 

Let  P(t)  be  a  solution  of  the  coupled  Riccati  equations  (6.2.1)-(6.2.2), 

P{t)  =  -A'P{t)  -  P(t)A  +  P{t)BxB[P(t)  -  C[CU  t  ±  %T  (6.3.46) 
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P{iT+ )  =  P(iT)  +  (B'2PiiT))'(j2I  -  B'2P(iT)B2)-1B,2P(iT).(6.3A7) 


It  can  be  rewritten  as  a  Lyapunov- type  equation, 

P{t)  =  —Ac(t)'P(t)  -  P(t)Ac(t)  -  H(t)'H(t),  t  ±  iT  (6.3.48) 
P(iT+)  =  P(iT)  +  Qtmp(iT).  (6.3.49) 

where  Ac(t),H(t)  and  QtmpiiT)  are  given  in  (6.2.17).  Therefore,  we  can  obtain 
the  following  results  from  [69,  71].  If  P(-)  is  a  symmetric  periodic  solution  of  the 
Riccati  equation  (6.3.46)-(6.3.47),  then  P(-)  is  also  a  solution  of  the  Lyapunov 
equation  (6.3.48)-(6.3.49).  The  structural  properties  of  the  pair  ( H(t),Ac(t )) 
can  be  related  to  the  ones  of  the  pair  (Ci,  A)  by  means  of  the  following  Lemma 
[71,  72]. 

Lemma  6.3.6  The  pair  (C\,  A)  is  observable  if  and  only  if,  for  any  T— periodic 
matrix  P(-),  the  pair  ( H(t),Ac(t ))  is  observable.  □ 

The  solution  of  equations  (6.3.48)-(6.3.49)  can  be  also  given  by  the  celebrated 
formula  [65].  Let  P0  be  the  initial  condition  of  equation  (6.3.48)  at  time  t0>  the 
solution  of  the  equation  (6.3.48)  is  given  by 

P(t)  =  &Ac(t0,t)Po$Ac(to,t)  -  f  &Ac(T,t)H(r)'H(T)&Ac(T,t)dT.  (6.3.50) 

Without  loss  of  generality,  we  assume  t0  =  0,  then  P(T )  is  given  by 

P(r)=4.'A,(r,0)-1P„^(r,0)-1-/T^(r.r)H(r)'if(T)®X(T,T)<iT.(6.3.51) 

J  0 

Since  we  are  looking  for  the  periodic  solutions  of  equations  (6.3.48)-(6.3.49),  the 
periodic  generator  Pq  must  satisfy  the  following  algebraic  Lyapunov  equation 

P0  =  &Ac(T,0)P0$ac(T,0)  +  W(T)  (6.3.52) 
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expressed  as 


N 

x(t)  —  Ax(t)  +  Biu(t)  +  B2  ^2  vk${t  ~  kT0) 

k= 0 

z(t)  =  Cx(t)  +  Du(t)  (5.3.52) 

where  S(t)  is  the  standard  Dirac  delta  distribution.  The  control  objective  is 
given  by 

II*WIIl,  <  72|M£.  (5.3.53) 

Remark  5.3.4  The  linear  model  has  a  physical  interpretation.  Since  each  im¬ 
pact  duration  is  very  short,  we  can  ignore  the  finite  duration  and  treat  the  impact 
as  an  impulsive  event.  This  linear  model  has  been  used  in  various  impact  prob¬ 
lems. 

We  should  point  out  that  the  control  problem  (5.3.52)-(5.3.53)  is  also  different 
from  the  standard  linear  ones  discussed  in  [20,  22].  The  disturbance  v  €  l2  is  a 
discrete-time  signal,  while  the  controlled  output  z(t)  are  continuous  signals.  The 
standard  approaches  [20,  22]  cannot  be  applied  to  this  problem  directly.  We  will 
analyze  and  solve  this  special  linear  control  problem  in  the  next  chapter. 
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Chapter  6 


#«,  CONTROL  FOR  IMPULSIVE  DISTURBANCES 


In  chapter  5,  we  studied  the  nonlinear  control  of  impact  forces  when  the  impact 
forces  are  modeled  as  finite  duration  events.  We  have  shown  that,  in  certain 
cases,  by  taking  advantage  of  the  fact  that  the  impact  duration  is  very  short, 
we  can  approximate  the  system  dynamics  undergoing  impact  by  low  order  terms 
of  the  Taylor  polynomials.  For  the  disturbance  attenuation  problem  the  leading 
term  provides  a  very  good  approximation.  A  linear  model  can  be  obtained  if 
we  treat  the  finite  event  by  an  impulsive  event.  Thus,  the  nonlinear  control 
problem  can  be  solved  via  a  linear  one.  In  this  chapter,  we  analyze  and  study 
the  disturbance  attenuation  problem  for  the  corresponding  linear  system.  This 
control  problem  has  some  special  features;  the  disturbance  input  v  €  I2  is  a 
discrete-time  signal,  while  the  controlled  output  z  is  a  continuous  signal.  It 
cannot  be  solved  by  the  existing  H0 0  approaches  [19,  20,  24,  21,  22,  23].  We  are 
motivated  by  the  approaches  [20,  24]  for  a  standard  H <*,  control  problem,  where  a 
state-space  solution  of  control  is  closely  related  to  the  solutions  of  the  Riccati 
equations.  If  full  state-feedback  is  available,  we  knew  that  a  controller  exists  if 
and  only  if  the  unique  solution  of  the  associated  algebraic  Riccati  equation  is 
positive  definite.  In  addition,  a  formula  for  the  state-feedback  gain  matrix  was 
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given  in  terms  of  the  solution  of  the  Riccati  equation.  In  this  control  problem, 
instead  of  the  one  algebraic  Riccati  equation  involved  in  problems  [20,  24],  we 
will  have  two  coupled  Riccati  equations  in  order  to  account  for  the  impulses 
in  the  state-feedback  case.  Hamiltonian  theory  is  used  to  analyze  the  coupled 
Riccati  equations,  and  necessary  and  sufficient  conditions  are  obtained  for  the 
existence  of  a  unique  solution.  Because  the  state  variables  may  experience  jumps 
due  to  the  effects  of  impulsive  forces,  an  extended  Lyapunov  lemma  is  derived 
to  prove  the  stability  of  the  closed-loop  system. 

6.1  Formulation  of  the  Control  Problem 

We  consider  a  finite-dimensional  linear  system 

z  Ln  L12  v 

y  I/2i  L22  u 

where  z(t),y(t)  are  the  controlled  and  measured  output  vectors,  respectively. 
They  are  piecewise-continuous  signals.  There  are  two  kinds  of  input  signals, 
u,  the  control  vector  which  is  assumed  to  be  piecewise-continuous,  and,  v,  the 
impulsive  disturbance  vector  which  is  assumed  to  be  in  the  Z2  space.  = 

1.2  are  linear  operators  mapping  from  {v,  u}  to  {z,y}.  This  is  a  hybrid  system 
because  it  contains  both  continuous  time  and  discrete  time  components.  In  view 
of  the  results  obtained  in  chapter  5  we  assume  that  the  system  (6.1.1)  admits  a 
state  space  realization  of  the  form 

OO 

x(t)  =  Ax(t)  -I-  Biu(t )  +  ^2  B2v(i)6(t  —  iT), 

i= 0 

z(t)  =  Cix{t )  4-  Diu(t), 

y(t)  =  C2x(t).  (6.1.2) 


(6.1.1) 
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It  is  convenient  to  use  the  state  space  equations  (6.1.2)  in  the  form 


x(t)  =  Ax{t)  4-  Biu(t),  t  ^  iT, 

x(t+)  —  x(t)  +  B2v(i),  t  =  iT,  i  =  0,1,2, . , 

z(t )  =  C\x(t )  +  Diu(t), 

y(t)  =  C2x(t).  (6.1.3) 


where  x(t+ )  denote  the  values  of  state  variables  immediately  after  a  jump.  The 
state  variables  x(t)  are  right  continuous  and  may  be  left  discontinuous  due  to 
the  possible  jump. 

The  control  problem  here  is  to  design  a  stabilizing  controller  to  attenuate  the 
effects  of  the  impulsive  disturbance  v  on  the  controlled  output  z.  Let  us  define 
K  as  the  set  of  all  causal,  finite-dimensional  linear  stabilizing  controllers. 

We  will  now  introduce  a  minimax  performance  measure  which  is  motivated 
by  [27].  For  k  G  K  define  a  performance  measure 


m  = 


sup  ( 

v€h,v^0 


Ml2) 

MV 


(6.1.4) 


J(k)  can  also  be  viewed  as  the  induced  operator  norm  from  l2  —>  L2. 

We  want  to  find  a  controller  k  G  K  to  minimize  the  worst  case  performance 
measure  J(k).  Specifically,  we  solve  the  following  (sub)optimal  problem.  Given 
7  >  0,  find  a  k  €  K  such  that  following  inequality  holds: 

J{k )  <  7 2  (6.1.5) 


Equivalently,  we  require  that 

<  0,  (6.1.6) 

for  all  possible  v  G  l2,  under  the  constraints  of  the  system  equations  (6.1.3).  If 
such  a  controller  k  exists,  we  call  it  a  7— level  disturbance  attenuation  controller. 
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6.2  Main  Results 


We  start  with  a  standard  H control  problem  treated  in  [20,  24]  where  we 
consider  a  continuous  LTI  system.  It  is  well  known  that  a  state-space  solution 
to  the  Jfoo  control  problem  is  closely  related  to  Riccati  equations.  If  full  state- 
feedback  is  available,  a  controller  exists  if  and  only  if  the  unique  solution  of  the 
associated  algebraic  Riccati  equation  is  positive  definite.  In  addition,  a  formula 
for  the  state-feedback  gain  matrix  was  given  in  terms  of  the  solution  of  the 
Riccati  equation.  Similar  results  can  be  obtained  for  discrete-time  systems  [22] 
and  time- varying  systems  [61]. 

It  should  be  noted  that  the  control  problem  defined  in  (6.1.6)  is  different  from 
the  standard  H ^  problems  treated  in  [20,  22,  24,  61].  Since  the  system  (6.1.3)  is 
a  hybrid  system  which  contains  both  continuous  and  discrete  components.  The 
control  problem  cannot  be  solved  by  the  formulas  obtained  in  [20,  24,  22,  61]. 
However,  the  essential  ideas  can  be  carried  over  to  analyze  the  control  problem 
(6.1.6).  For  state-feedback  control,  instead  of  the  one  algebraic  Riccati  equation 
involved  in  problems  [20,  24],  we  will  have  two  coupled  Riccati  equations  given 
by 

K(t)  =  -A! K{t)  -  K(t)A  +  K{t)BiB[K{t)  -  C[CU  t  ±  iT  (6.2.1) 
K(iT+)  =  K(iT)  +  {B'2K(iT))'(kl  ~  B'2K{iT)B2)-1B'2K(iT)  (6.2.2) 

where  7  >  0  is  a  given  real  number,  the  n  x  n  matrix- valued  function  K{t)  is 
right  continuous  and  may  be  left  discontinuous,  K(iT+)  represents  the  value  of 
K{t)  after  the  ith  jump,  while  K{iT)  represents  the  value  of  K(t)  just  before  the 
ith  jump,  i.e.,  K(iT )  :=  lim£>0ie->o  K(iT  —  e)  assuming  that  the  limit  exists. 

The  main  results  of  this  chapter  are  stated  as  Theorems  6.2.1-  6.2.3.  We 
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assume  that  all  state  variables  are  measured,  and  consider  the  infinite  horizon 
problem  with  zero  initial  state.  First  the  following  assumptions  are  made. 

i)  (A,  B\)  is  controllable,  {Ci,A)  is  observable. 

ii)  D'fiC,  Di]  =  [0  /]. 

Assumption  i)  is  standard  in  the  quadratic  regulation  of  a  linear  system. 
It  can  be  relaxed  by  the  assumption  that  ( A,B\ )  is  stabilizable  and  (C\,  A)  is 
detectable.  Assumption  ii)  is  made  here  just  for  the  sake  of  simplicity.  Relaxing 
this  assumption  will  only  complicate  the  formulas,  but  an  analysis  can  be  carried 
out  along  lines  similar  to  what  appears  below.  We  will  state  the  main  theorems 
first  and  defer  the  proofs  to  the  following  sections. 

Theorem  6.2.1  Consider  the  hybrid  system  described  by  (6.1.3).  Let  7  >  0  be 
given.  Suppose  that  the  assumptions  i )  and  ii)  hold.  Then  there  exists  a  con¬ 
troller  k  €  K  such  that  J(k )  <  72  if  there  exists  a  unique  stabilizing  positive 
definite  periodic  solution  P(t)  of  the  coupled  Riccati  equations  (6.2.1)-(6.2.2). 
Moreover,  if  this  condition  is  satisfied,  one  such  stabilizing  state-feedback  con¬ 
troller  is  given  by 

u(t )  =  -B[P(t)x(t),  Vt  >  0.  (6.2.3) 

□ 

Remark  6.2.1  As  7  — »  00,  the  coupled  Riccati  equations  (6.2.1)-(6.2.2)  degen¬ 
erate  into  the  continuous-time  Riccati  equation 

k{t)  =  —A'K(t)  —  K(t)A  +  K(t)BiB[K(t)  —  C[C\.  (6.2.4) 

This  Riccati  equation  will  yield  a  unique  positive  definite  constant  matrix  solu¬ 
tion  under  assumption  i ).  It  is  well  known  that  this  unique  solution  internally 


92 


stabilizes  the  associated  closed-loop  system  from  LQG  theory.  Thus,  the  H0 0 
problem  degenerates  into  a  LQG  problem  as  7  — >  00. 

Via  the  Hamiltonian  theory,  the  solution  of  a  Riccati  equation  can  be  obtained 
by  solving  a  suitable  set  of  linear  differential  equations.  The  Hamiltonian  matrix 
associated  with  the  continuous  time  Riccati  equation  (6.2.1)  is  given  by 

A  —B1B[ 

-C'XCX  -A' 

As  usual  let  us  consider  a  2 nth  order  differential  equation 

X(t)  =  HX(t),  t^iT.  (6.2.6) 

The  state  transition  matrix  associated  with  H  is 


(6.2.5) 


X(t)  =  $#(£, t)X(t),  t>r ,  and  t,r^iT,  (6.2.7) 

where  $#(£,  r)  has  the  following  properties, 

d^,T)  =  r),  $H(r,  r)  =  I.  (6.2.8) 

The  symplectic  matrix  associated  with  the  difference  Riccati  equation  (6.2.2)  is 
given  by 

I  -r2B2B'2 

F  = 

0  / 

Let  us  define  a  2 nth  order  difference  (jump)  equation 

X(t+)  =  FX(t),  t  =  iT.  (6.2.10) 

The  combined  equations  (6.2.6)  and  (6.2.10)  are  a  hybrid  system.  For  r  =  iT,i  — 
0, 1,  **•,***,  and  r  <  t  <  (i  +  1)T,  we  have, 

X(t)  =  $H(t,T)FX(r), 

=  &(t,T)X(r).  (6.2.11) 


(6.2.9) 
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We  will  show  later  that  this  state  transition  matrix  $(£,  r)  of  the  hybrid  system 
displays  the  same  properties  as  for  a  periodic  system. 

Now,  let  P(t)  be  a  periodic  solution  of  equations  (6.2.1)-(6.2.2).  The  associ¬ 
ated  closed-loop  system  is  given  by 

x(t)  =  [A  —  BiB[P(t)\x(t).  (6.2.12) 

Definition  6.2.1  A  periodic  solution  P(t)  of  (6. 2.1)-(6. 2.2)  is  called  a  stabiliz¬ 
ing  solution  if  the  closed-loop  system  (6.2.12)  is  asymptotically  stable. 

Let  us  consider  the  following  hybrid  system  first: 


x{t)  =  [A-BxB[P(t)]x(t), 

=  Ac(t),  t  iT , 

(6.2.13) 

x(t+)  =  [I-y-2B2B'2P(t)]x(t), 

=  Fc(t),  t  =  iT. 

(6.2.14) 

We  will  show  that  if  (6.2.12)  is  asymptotically  stable,  then  the  hybrid  system 
asymptotically  stable  also. 

Theorem  6.2.2  Let  7  >  0  be  given.  Suppose  that  the  assumptions  i)  and  ii) 
hold.  Then  a  necessary  and  sufficient  condition  for  the  existence  of  a  unique 
positive  definite  periodic  solution  P(t )  to  equations  (6.2.1)  -  (6.2.2),  such  that 
the  hybrid  system  (6.2.13)-(6.2.14)  is  asymptotically  stable,  is  that  no  eigenvalue 
o/$(T,  0)  lies  on  the  unit  circle.  □ 

The  next  theorem  gives  necessary  and  sufficient  conditions  that  the  unique  pos¬ 
itive  definite  periodic  solution  stabilizes  the  associated  closed-loop  system. 
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Theorem  6.2.3  Suppose  that  the  assumptions  i )  and  ii )  hold.  Let  P(t )  be  the 
unique  positive  definite  periodic  solution  of  (6.2.1)  and  (6.2.2).  Then  P(t )  is  a 
stabilizing  solution  if  and  only  if  the  inequality 

WJX)  -  «^(T,0)Qtmp(T)$/1„(T,0)  >  0.  (6.2.15) 

holds,  where 

WC(T)  =  [T  $Ac(T,0)'H(T)'H(T)<t>Ac(T,0)dT  (6.2.16) 

J  0 

and  $Ac(t,T)  is  state  transition  matrix  of  Ac(t),  and  where  the  matrices  Ac(t), 
H(t)  and  Qtmp(T)  are  defined  by 

Ac(t )  =  A  —  BiB[P(t), 

H(t)’H(t)  =  C[C\  +  P(t)BiB[P(t), 

QtmV(T)  =  (B‘2P(T))'('I2I-B'2P(T)B2)-iB^P[T).  (6.2.17) 

□ 


6.3  Background  and  Technical  Lemmas 

In  order  to  prove  Theorems  (6. 2.1)- (6. 2. 3),  we  need  to  develop  theory  and  tech¬ 
nical  machinery  to  analyze  the  coupled  Riccati  equations  (6.2.1)-(6.2.2).  Specif¬ 
ically,  1)  we  will  give  conditions  for  existence  of  periodic  solutions  of  equations 
(6.2.1)-(6.2.2),  2)  If  such  conditions  are  satisfied,  we  will  parameterize  all  stabi¬ 
lizing  periodic  solutions.  The  technical  machinery  developed  here  is  based  on  the 
standard  tools  for  the  analysis  of  periodic  systems  [62,  63,  64].  We  will  show  that 
equations  (6.2.1)-(6.2.2)  display  same  properties  similar  to  a  standard  periodic 
system,  the  main  difference  being  that  we  need  to  take  care  of  the  jumps  in  the 
state  which  actually  result  in  the  periodic  solutions. 
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6.3.1  Periodic  system  and  stability  analysis 

Consider  the  following  linear  periodic  system 

x(t)  =  A(t)x(t )  +  B(t)u(t ) 

y(t)  =  C(t)x(t)  (6.3.1) 

where  x(t)  G  Rn,y(t)  G  Rm,u(t)  G  RP  are  state,  output  and  input  vector  vari¬ 
ables  respectively,  A(-),  £?(•),  C(-)  are  assumed  to  be  periodic  matrices  of  period 
T,  i.e, 


A(t  +  T)  —  A(t),  B(t  +  T)  =  B(t),  C(t  +  T)  =  C(t ),  W  >  0.  (6.3.2) 

The  state  transition  matrix  of  system  (6.3.1)  is  denoted  by  ^(t,t0).  The  matrix 
ty(t  +  T,t)  is  called  the  monodromy  matrix  at  time  t.  Since  4>(f  +  T,  t)  = 
^>(t,  t)^(t  +  T,  r)^(t,  r)_1,  the  eigenvalues  of  ^( t  +  T,  t)  are  independent  of  t. 
they  are  also  called  the  characteristic  multipliers  [65,  66,  67]. 

Remark  6.3.1  [67,  68]  The  system  (6.3.1)  is  asymptotically  stable  if  and  only 
if  the  eigenvalues  of^(T,  0)  are  all  inside  the  unit  circle. 

Definition  6.3.1  [62]  An  eigenvalue  A  ofty(T,  0)  is  said  to  be  ( C(t),A(t ))  un¬ 
observable  if  there  exists  an  n-dimensional  vector  £  ^  0  such  that 

*(T,  0)£  =  A£  (6.3.3) 

C't^O)^  =  0,  Vt  G  [0,  T].  (6.3.4) 

Remark  6.3.2  The  pair  ( C(t),A(t ))  is  said  to  be  detectable  if  and  only  if  all 
eigenvalues  A  o/^(T,  0)  such  that  |A|  >  1  are  ( C(t),A(t ))  observable. 


96 


at  the  sampling  instances  are  given  by 


x(k  +  l)  =  Fx(k)  +  Gxu(k)  +  G2v(k) 

y(k)  =  C2x(k)  (6.5.2) 

where  F  =  eAT,  G\  =  /0T  eA^T~T>>  B\dr .  G2  ~  eATB2.  The  controlled  output  is 
continuous  and  is  given  by 

z(kT  +  t)  =  C\eAtx{kT)  +  /*  BxdTu(kT)  +  CxeAtB2v{kT ) 

Jo 

Vi  6  [kT,  kT  +  T}.  (6.5.3) 

It  is  well-known  if  (A,  Bx )  is  controllable  and  (C{,  A)  is  observable,  for  a  generic 
choice  of  the  sampling  period  T  so  are  the  pairs  ( F,Gi )  and  ( CX,F )  [74].  We 
assume  that  the  sampling  period  T  satisfies  generically  condition  in  the  rest  of 
this  chapter. 

Let  K  be  the  set  of  all  discrete-time  LTI  stabilizing  controllers.  The  control 
problem  we  address  here  is  that  for  given  7  >  0,  find  a  k  €  k  such  that  the 
following  inequality  holds 

k(K)  =  IMI12-72|«  <  0.  (6.5.4) 

for  all  v  £  l2,  and  zero  initial  conditions.  It  is  known  that  a  sampled-data 
control  system  in  general  is  time-varying  [76].  H0 0  control  design  problem  for 
sampled-data  control  systems  has  been  addressed  by  many  researchers  recently  ( 
see  [23,  77,  75]  for  details).  Our  interest  here  is  to  adapt  the  theory  to  the  case  of 
impulsive  disturbances  attenuation.  We  will  show  that  by  using  a  sampled-data 
controller,  the  control  problem  can  be  converted  into  a  LTI  discrete-time 
control  problem. 
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The  term  \\z\\\2  during  [kT,  ( k  +  1)T]  can  be  rewritten  as 
fT 

/  z'(kT  +  t)z{kT  +  t)dt  =  x' {k)M\x{k)  +  2u' (k)P[x(k)  +  u' (k)M2u(k)  + 

J  0 

2v'(k)P^x(k)-hj\k)M3v(k)+2v'{k)P^u(k).  (6.5.5) 

where  Mj  are  weighting  matrix  associated  with  state  x(k),  control  u(k)  and 
disturbance  v(k),  and  Pt  are  cross-terms  respectively.  They  are  given  by, 

Mx  =  F  eA'tC[CleMdt,  P1=fT  te-W^drdt 
Jo  Jo  Jo 

M2=  [  [' B'^'^C'.C^-^B.drdt,  P2=fT  eA'tC[C1eMB2dt 
Jo  Jo  Jo 

M3=  [T  By^C^B^t,  P3=/r  f  By  tC,xCyt~^B2drdt.  (6.5.6) 
Jo  Jo  Jo 

Hence  the  square  of  the  controlled  output’s  norm  is  in  turn  given  by 

roo  °°  p(k+l)T 

INI La  =  L  Z'(t)z(t)dt  =  E  /  Z'{t)z(t)dt 

Jo  k=0Jk 

oo  fT 

=  z{kT  +  t)'z{kT  +  t)dt 

k= o  Jo 
00 

=  ^2  x'(k)Mix(k)  +  2u'(k)P[x(k)  +  v!(k)M2u{k)  + 

fc=0 

2v\k)P!2x{k)  +  v\k)Mzv{k)  +  2  v'{k)Pyk).  (6.5.7) 

The  sampled-data  control  problem  is  now  converted  into  the  following  discrete¬ 
time  Poo  control  problem,  where  the  system  is  given  by, 

x(k  +  1)  =  Fx(k)  +  Giu(k)  +  G2v(k) 

y(k)  =  C2x(k),  k  =  0, 1,2,  •••  (6.5.8) 

and  the  control  objective  is  to  find  a  k  £  K  such  that 

J7(P)  =  |N||2  -  T2|NII?2  <  0  (6.5.9) 


119 


where  ||z|||2  is  given  by  (6.5.7). 

Among  different  state-space  approaches  to  this  class  of  H ^  control  problems, 
the  one  that  uses  the  framework  of  dynamic  game  theory  seems  to  be  the  most 
natural.  This  is  so  because  the  original  7/oo-optimal  control  problem  (in  its 
equivalent  state-space  formation)  is  in  fact  a  minimax  optimization  problem, 
and  hence  a  zero-sum  game  [27],  where  the  controller  can  be  viewed  as  the 
minimizing  player  and  disturbance  as  the  maximizing  player.  We  will  use  this 
approach  to  derive  the  (sub)optimal  solution  to  the  problem  (6.5.8)-(6.5.9). 

6.5.2  Full  state-feedback  control 

Consider  a  standard  discrete-time  system  given  by 

x{k  +  1)  =  Ax{k)  +  B\u{k )  +  B2v(k) 

y(k)  =  C2x(k).  (6.5.10) 

The  performance  measure  is  given  by 
00 

J7(K)  =  ^2  x'(k)Qx(k)  +  u'(k)u(k)  —  7 2v'(k)v(k)  <  0  (6.5.11) 

k= 0 

for  some  specified  7  >  0  and  Q  >  0,  and  v  G  l2- 

First  we  consider  the  case  when  the  full  state  variables  (C2  =  I )  are  available. 
We  make  the  following  standard  assumptions: 

Al)  The  pairs  (A,  B\)  and  (A,  B2 )  are  controllable. 

A2)  (A,  Q1/2)  is  observable. 

Let  us  define  an  algebraic  Riccati  equation, 

M  =  A'{M~l +BlB[-~i-2B2B'2)-lA  +  Q.  (Rl) 
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Theorem  6.5.1  [27]  Consider  discrete-time  system  described  by  (6.5.10),  let 
7  >  0  be  given.  Suppose  that  the  assumptions  Al )  -  A2 )  hold.  Then  there  exists 
a  controller  k  e  K  such  that  inequality  (6.5.11)  hold  if  and  only  if 

1)  (Rl)  admits  a  positive  definite  solution  M+. 

2)  7 2I  -  B'2M+B2  >  0. 

Moreover,  when  these  conditions  hold,  one  such  stabilizing  feedback  controller  is 
given  by 


u{k) 

=  -B[M+&-lAx{k), 

(6.5.12) 

v(k) 

=  r2B'2M+A~lAx(k ), 

A  : 

=  I  +  (BiB[  -7-2B2B')M+.  A:  =  1,2, 3, . 

• 

and  the  closed-loop  matrix  ( I  —  BiB[M+  A  is  Hurwitz.  □ 

Now  we  are  ready  to  present  our  main  results.  Let  us  first  make  the  following 
assumptions: 

Al)  7 2 1  -  M3  >  0, 

A2 )  The  pairs  (A,B{),  ( A,B2 )  are  controllable  and  (Ci,A)  is  observable, 

A3)  Bi  is  full  rank, 

A4 )  Q  >  0  and  is  defined  by 

Q  =  Mi  +  P2MflP2  ~  (Pi  +  PsM^P^'^iPi  +  PzM^Pf)  (6.5.13) 
where  the  M2  and  M3  are  given  by 

M3  =  72/  —  M3,  M2  =  M2  +  T3M3  P%.  (6.5.14) 

and  Mi  and  Pi  are  defined  by  (6.5.6). 
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Theorem  6.5.2  Consider  the  sampled-data  system  described  by  (6.5.1).  Let 
7  >  0  be  given.  Then  there  exists  a  controller  k  €  K  such  that  (6.5.9)  hold  if 
and  only  if  following  conditions  hold: 

1’)  The  Riccati  equation  (Rl)  associated  with  (A,  Bi,  B2,Q)  has  a  positive  defi¬ 
nite  solution  M+, 

2’)  7 2 1  -  B2M+B2  >  0. 

When  these  conditions  hold,  one  such  stabilizing  controller  is  given  by 

u(k)  =  -Bi  M+A~xAx(k),  (6.5.15) 

v(k)  =  —'y~2B2M+A~1Ax(k), 

A:  =  /+(B15/-7-2B2B2')M+.  k  =  1,2,3, . 

where  the  A,BX,B2  are  defined  by, 

A  =  F  +  G1M3“1P2-  hGlMf1Pi)Mfl(P[  +  PzMflP!fi) 

Bx  =  {G1  +  'yG2Mf1P^)M2-1/2 

B2  =  7  G2Mf1'2  (6.5.16) 

Proof:  In  order  to  apply  theorem  (6.5.1)  to  this  problem,  we  first  need  to  con¬ 
vert  the  problem  (6.5.8)-(6.5.9)  into  the  standard  problem  described  by  (6.5.10)- 
(6.5.11),  then  we  show  that  the  assumptions  of  theorem  (6.5.1)  are  satisfied. 
Since  the  pairs  (A,  B\)and{A,  B2)  are  controllable  and  ( Ci,A )  is  observable  by 
assumption  A2),  hence  the  discrete  system  (F,G i),  (F,  G2)  are  also  controllable 
and  (Ci,  F)  is  observable  by  the  assumption  of  that  generic  condition  of  the  sam¬ 
pling  period  T  is  satisfied. 

Next,  we  carry  out  linear  transformations.  By  doing  so,  we  first  define  a  new 


122 


variable  v(k)  as 

v(k)  = -M3l^2[v(k)  -  M3~1(P2x(k)  +  P^u(k))],  k  =  0, 1,2,  •••  (6.5.17) 
7 

The  M3  is  invertible  by  assumption  AT).  The  system  (6.5.8)  and  the  performance 
measure  (6.5.9)  become 

x(k  +  1)  =  (A  +  G2M3~lP!f)x{k)  +  (Gx  +  G2MflP!f)u(k)  +  ~/G2M3~1/2v(k) 

OO 

Jy{K)  =  Y,x\k)(Mr+P*MflP2)x(k)  +  2u{k),{Pl  +  PzMflP!1)xtk)  + 

k=0 

u'(k)(M2  +  P3MflP'z)U(k)  -  7  2v(k)'v(k) 

Since  (C\,A)  is  observable,  and  B\  is  full  rank  by  assumption  A3),  it  is  easy  to 
show  that  the  matrix  M2  is  positive  definite.  The  M3~l  >  0  follows  from  M3  >  0. 
This  implies  that  P3M3~1P^  >  0.  Hence,  M2  =  M2  +  P3M3  1P^  >  0  is  positive 
definite.  Define  a  new  variable  u(k)  as 

u(k)  =  M21/2[u(k)  +  M2~\P[  +  P3M3~lP!ff\x{k)  (6.5.18) 

andQ  =  M1  +  P2M3~1P^-(P{  +  P3Mf1P!l),M2~1{P{  +  P3M3~1P^).  After  some 
simplifications,  we  finally  arrive  at 

x(k  +  1)  =  Ax(k)  +  Biv[k)  +  B2u(k)  (6.5.19) 

OO 

J-y(K)  =  ^2  x(k)'Qx(k)  +  u(k)'u(k)  —  7 2v(k)'v(k)  (6.5.20) 

fc= 0 

Since  the  linear  transformations  do  not  affect  the  controllability  and  observabil¬ 
ity,  it  follows  that  (A,  Bi)and(A,  B2)  are  controllable.  Q  is  positive  semi-definite 
by  its  structure  and  assumption  AA).  Since  ( C\,A )  is  observable,  it  is  straight¬ 
forward  to  show  that  (A,  Q1^2)  is  also  observable.  Hence,  the  assumptions  AT) 
and  A2)  of  theorem  (6.5.1)  are  satisfied,  and  the  proof  is  completed  by  invoking 
the  theorem  (6.5.1)  to  the  problem  (6.5.19)-(6.5.20).  □ 
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The  feedback  controller  is  expressed  by  in  the  new  variables  u  and  v.  By  changing 
back  the  variables,  we  get 

u(k)  =  -M^i^M+A-'A-M^iPi  +  PsM^P^xik),  (6.5.21) 
v(k)  =  -M^l/2[^-lB'2M+/\-lA  -  Mf)1/2P^]x{k)  +  M^P^u{k). 

In  most  practical  control  applications,  only  partial  state  variables  are  available, 
we  will  derive  the  solution  for  output-feedback  problem  in  next  section. 

6.5.3  Output-feedback  control 

When  the  full  state  variables  are  not  available,  we  need  to  construct  a  observer- 
based  feedback  controller.  The  following  theorem  gives  us  the  standard  results. 
We  first  assume  the  following  assumptions  hold. 

Al)  The  pairs(A,Pi),  (A,  B2)  are  controllable  and  (C2,  A)  is  observable. 

A2)  Q  is  positive  semi-definite. 

Let  us  define  two  algebraic  Riccati  equations  as  follows, 

M  =  A'iM-1 +B1B[-j-2B2B'2)-1A  +  Q  (PI) 

N  =  A(N~X  +  C2C2  —  ry~2Q)~lA!  +  B2B'2  (R2) 

Theorem  6.5.3  [27]  Consider  discrete-time  system  described  by  (6.5.10),  let 
7  >  0  be  given.  Suppose  that  the  assumptions  Al)  —  A2)  hold.  Then  there  exists 
a  controller  k  G  K  such  that  inequality  (6.5.11)  hold  if  and  only  if 

1)  (Rl)  and  (R2)  have  positive  definite  solutions  M+  and  N+. 

2)  p(M+N+ )  <  72. 

Moreover,  when  these  conditions  hold,  one  such  stabilizing  feedback  controller  is 
given  by 

x(k  +  1)  =  Ax(k)  +  A(A+-1  +  C'2C2  -  'y~2Q)~1[y~2Qx(k) 
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+C'2(y(k)  -  C2x(k))\  +  Biu(k)  (6.5.22) 

+  B&  -  r2B2B'2)-1A(I  -  yiN+M+)-1x(k), 


u(k)  = 

where  the  x  is  an  observer.  □ 

Now  we  solve  the  problem  (6.5.8)-(6.5.9)  by  applying  theorem  (6.5.3).  Let  us 
first  make  the  following  assumptions: 

AT)  7 2 1  -  M3  >  0, 

A2)  The  pairs  (A,  BT)  is  controllable  and  {C\,A)  is  observable,  so  are  (A,B2) 
and  (C2,A)  respectively, 

Ad)  Bi  is  full  rank, 

A4)  Q  >  0  and  is  defined  by 

Q  =  MX  +  P2M3_1P2  -  [P[  +  PzMflP2)'MAl(P[  +  PzMflP2)  (6.5.23) 
where  the  M2  and  M3  are  given  by 

M3  =  72/  —  M3,  M2  =  M2  +  P3M3  P3.  (6.5.24) 

and  Mi  and  Pi  are  defined  by  (6.5.6). 

Theorem  6.5.4  Consider  the  sampled- data  system  described  by  (6.5.1).  Let 
7  >  0  be  given.  Then  there  exists  a  controller  k  G  K  such  that  (6.5.11)  hold  if 
and  only  if  following  conditions  hold: 

1’)  The  two  Riccati  equations  (Rl)  and  (R2)  associated  with  (A,  Bi,  B2,C2,Q) 
have  positive  definite  solutions  M  and  N. 

2’)  p(MN)  <  72. 

When  these  conditions  hold,  one  such  stabilizing  controller  is  given  by 
x(k  +  1)  =  Ax(k)  +  A(N-'  +C2C2-T'2Q)~1b~2Qx{k) 
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+C2(y(k)  -  C2x(k))]  +  Biu(k)  (6.5.25) 

-BiiM-1  +  BxBi  -  7 -2B2B2')-1A(I  -  72iVM)-1x( k) 


u(k)  = 

where  the  (A,  B\,  B2,  C2)  are  defined  by, 

A  =  F^G2MflP!2-{1G2MflPl)M2-\P[  +  PzMflP^) 

Bl  =  (G1  +  7G2M3_1P3)^2_1/2 

b2  =  7g2m3“1/2 

C2  =  C2.  (6.5.26) 

Proof;  The  proof  is  analogous  to  the  proof  of  theorem  (6.5.2)  and  is  omitted 
here.  □ 

The  output-feedback  controller  (6.5.25)  is  expressed  by  in  the  new  variable  u(k), 

and  the  actual  control  u(k)  can  be  obtained  by  changing  back  the  variable 
(6.5.18) 

£(*  +  1)  =  Ax{k)  +  A(N~1  +  C2C2  —  ry~2Q)~1['y~2Qx(k)  +  C2  (y(k) 

-C2x{k))}  +  BiM2~1/2(P[  +  P3M3_1P^)x(fc)  +  BxM2/2u{k) 
u(k )  =  -M2~V2Bx{M~l  +  BxBi  -  'y~2B2B2)-1A{I  -  72iVM)-1x( k) 
-M2~\p[  +  PzMflP£)x(k).  (6.5.27) 

To  recap,  let  us  list  the  steps  in  the  design  procedures  for  the  output  feedback 
problem: 

step  1:  Calculate  the  matrices  M;  and  Pi  according  to  (6.5.6), 

step  2:  Give  7  >  0,  check  if  the  conditions  Al)  and  AA)  are  satisfied,  if  not, 

increase  7  until  these  conditions  are  satisfied, 

step  3:  Calculate  matrices  A,  B\,  B2,  C2  according  to  (6.5.26), 
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step  4:  Solve  two  algebraic  Riccati  equations  (Rl),  ( R2 )  for  M,  N  associated  with 
(A,  Ui,  B2 ,  C2,  Q ),  check  if  conditions  1'  and  2'  are  satisfied,  if  not,  increase  7, 
go  back  to  step  3. 

The  state-feedback  is  a  special  case,  the  design  procedure  is  similar  to  output 
feedback. 


6.5.4  Example 


Our  motivation  for  developing  the  Hqo  optimal  control  for  impulsive  disturbances 
is  from  various  impact  control  problems.  Let  us  consider  the  following  impact 
control  example  shown  in  Figure  6.2.  The  system  represents  a  flexible  body  with 
mass  m2,  and  the  k,  c  represents  its  stiffness  and  damping  coefficient  respectively, 
and  u  is  the  control  input.  The  system  is  subjected  to  a  series  of  impulsive 
disturbances  (impact  forces).  When  the  impact  occurs,  a  significant  portion 
of  the  kinetic  energy  of  the  impactor  will  be  transferred  to  the  flexible  body, 
causing  vibrations.  Hence,  a  reasonable  control  objective  is  to  minimize  the 
energy  transferred  to  the  system  by  the  impact  forces.  This  control  problem 
can  be  viewed  as  a  disturbance  attenuation  problem.  Suppose  the  displacement 
of  the  mass  is  available.  The  state  space  representation  of  the  system  with 
impulsive  disturbances  is  given  by, 


x(t)  = 


0  1 
—k/m2  —c/m2 

0 

fc=°  I  b2/m2 


x(t)  + 


0 

1/  m2 


u{t) 


+  E 


v(k)S(t  -  kT) 


=  Ax{t)  +  B\u(t )  +  ^2  B2v(k)S(t  —  kT) 

k= 0 
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where  the  W  ( T )  is  given  by 

W(T)  =  WC(T)  -  4>'Ac(T, 0)Qtmp(T)<bAc (T, 0).  (6.3.53) 

Proposition  6.3.1  [70]  The  Lyapunov  equations  (6.3.48)-(6.3.49)  admit  a 
unique  positive  periodic  solution  if  and  only  if  the  algebraic  Lyapunov  equation 
(6.3.52)  admits  a  positive  definite  solution  P0.  □ 

Proposition  6.3.2  [68]  Suppose  that  the  Ac(t )  is  asymptotically  stable.  Then 
the  necessary  condition  that  algebraic  Lyapunov  equation  (6.3.52)  admits  a  pos¬ 
itive  definite  solution  is  that  W(T)  defined  in  (6.3.53)  is  positive  definite.  □ 

6.4  Proofs  of  Theorem  6.2.1-  6.2.3 

Now,  we  are  ready  to  prove  the  main  results  of  this  chapter.  First  we  prove  the 
Theorem  (6.2.2),  then  we  prove  the  Theorem  (6.2.3).  Finally,  we  consider  the 
Theorem  (6.2.1). 

6.4.1  Proof  of  theorem  6.2.2 

Recall  the  hybrid  system  defined  in  (6.2.13)-(6.2.14), 

x(t )  =  [A-  BxB[P{t))x{t ), 

=  Ac{t),  t^iT,  (6.4.1) 

x(t+)  =  [I  -  7 ~2B2B'2P(t)]x(t), 

=  Fc(t),  t  =  iT.  (6.4.2) 

Denote  the  state  transition  matrix  of  the  hybrid  system  (6.4.1)  -(6.4.2)  by 
$c(t,  r).  It  has  the  following  properties, 

-  Ac(t)®c(t,r),  (6.4.3) 
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$c(T+,T) 


(6.4.4) 


=  Fc(t). 

Theorem  6.4.1  Let  P(t)  in  equation (6. 3. 24)  be  substituted  into  equation 

(6. 4-1)  and  (6.4-2).  Also  let  X,i  =  1,2, . n  be  n  eigenvalues  o/$(T,  0)  used 

to  form  Ax  in  equation (6. 3. 12),  then  <&c(t,  0)  is  given  by 

*cM)  =  Y{t)Y{  0)-1  (6.4.5) 

and  the  eigenvalues  o/$c(T,  0)  are  \,i  =  1,2, . n. 

Proof;  Using  equation  (6.3.15),  we  have 

*cM)  =  mnor1 

=  [AY(t)  -  B^XWYiO)-1 

=  [A-B1B[P(t)}^c(t,0) 

=  Ac(t)<f>c(t,  0).  (6.4.6) 

$c(o+,o)  =  r(o+)y(o)-1 

=  [Y(0)-1-2B2B,2X(0)]Y(0)-1 

=  [I-r2B2B'2P{  0)] 

=  Fc(  0).  (6.4.7) 

From  equation  (6.3.25), 

Y(T)  =  y(0)Ax 

$c(t,  o)  =  ycr)y(o)-1  =  y^AxF^)-1.  (6.4.8) 

Thus,  eigenvalues  o/$c(T,  0)  are  \,i  =  1,2, . n.  □ 

The  following  Lemmas  relate  the  positive  definite  periodic  solution  of  the  coupled 
Riccati  equations  and  the  structure  properties  of  systems  (controllability  and 
observability) . 
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Lemma  6.4.1  Assume  (C\,A)  is  observable  and  no  eigenvalue  o/$(T,  0)  lies 
on  the  unit  circle.  Then  there  exists  a  positive  definite  periodic  solution  P(t)  to 
equations  (6.2.1)  -  (6.2.2)  such  that  the  hybrid  system  (6.4-l)-(6.4-2)  is  asymp¬ 
totically  stable  if  {A,  Bfi)  is  controllable. 

Proof;  Since  only  solutions  P(t)  of  equations  (6.2.1)  and  (6.2.2)  that  lead  to 
asymptotically  stable  hybrid  systems  are  under  consideration,  n  eigenvalues  in 

Ai  must  be  chosen  so  that  | A*|  <  1  fori  —  1,2, . ,  n  according  to  Theorem 

(6.4-1).  Such  a  choice  is  possible  under  the  current  hypothesis  and  Theorem 
(6.3.1).  From  Theorem  (6.3.3),  the  resulting  P(t)  clearly  satisfies  P(t)  >  0.  The 
non-existence  of  P(t)  >  0  then  arises  only  when  the  condition  that  |Y(£)|  0  for 

all  t  £  [0,  T]  is  violated.  In  the  sequel,  we  will  show  that  (A,  Bfi)- controllability 
implies  |F(£)|  ^  0  for  all  t  £  [0,T].  Suppose  to  the  contrary  that  for  some 
t  =  s,s  £  [0,  T],Y(t)  is  singular.  Let  b  ^  0  be  an  n-dimensional  vector  belonging 
to  the  null  space  ofY(s),  i,e. 

Y(s)b  =  0  (6.4.9) 

Also  for  k  —  0,1, . .  Let  a(k)  and  fi(k)  be  n-dimensional  vectors  such  that 

a(k) 

. m . 

Define  a  vector  u(k)  as 

u(k)  = 

From  equation  (6.3.15),  we  have 

u{k  +  l)  =  <3>(5  +  T,  S)u(k)  (6.4.12) 


Y(s)Ak1b 

X{s)Ak1b 


a(k) 

m 


(6.4.10) 


(6.4.11) 
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Let  us  define  a  scalar  W(k )  as 


W(k)  =  fi*(k)a(k)- /3*(k  +  l)a(k  +  l)  (6.4.13) 

Since  o(0)  =  Y(s)b  =  0,  and  a{k),fi(k)  — >  0  as  k  — >  oo, 

OO 

Y,W(k)  =  0  (6.4.14) 

k= 0 

On  the  other  hand,  from  equation  (6.3.39) 

W(k )  =  u*(k  +  l)Vu{k  +  l)-u*(k)Vu{k) 

=  u*(k)[$(s  +  T,  s)'y$(s  +  T,  s)  -  V]u(k) 

=  u*{k)M{s  +  T,s)u(k)  (6.4.15) 

Since  M(s  +  T,  s)  >0  from  equations  (6.3.43)-(6.3.45),  hence  yielding  W(k )  > 

0.  Equation  (6. 4-H)  then  implies  that  W(k)  =  0  for  all  k  =  0,1,2, . . 

Therefore,' M(s  +  T,  s)  is  not  positive  definite,  hence  (A,  B{)  is  not  controllable, 
a  contradiction.  □ 

Lemma  6.4.2  Assume  that  no  eigenvalue  of  ${T,  0)  lies  on  the  unit  circle  and 
(A,Bi)  is  controllable.  Then  there  exists  a  unique  positive  definite  periodic  so¬ 
lution  P(t )  to  equations  (6.2.1)-(6.2.2)  if{C\,A )  is  observable. 

Proof:  Under  the  assumptions,  the  existence  of  at  least  one  solution  P(t)  >  0  is 
guaranteed  by  Lemma  (6.4-1)-  Such  a  solution  P(t )  is  obviously  associated  with 
the  asymptotically  stable  hybrid  system  (6.4-1)  and  (6. 4-2).  In  the  following,  we 
shall  show  that  the  additional  condition  {C\,  A) -observability,  necessarily  implies 
the  asymptotic  stability  of  the  hybrid  system  (6-4- 1)  and  (6-4-2) .  The  uniqueness 
of  the  solution  can  be  established  by  this  way.  Assume  to  the  contrary  that  the 
hybrid  system  (6-4-1)  and  (6-4-2)  is  unstable.  It  then  follows  from  Theorem 
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(6. 4-1)  that  hi  in  equation  (6.3.12)  has  at  least  one  eigenvalue  A  o/4>(T,  0)  such 
that  |A|  >  1.  Since  the  order  of  eigenvalues  in  Ai  is  irrelevant  to  solution,  it  can 
be  assumed  that  Xx  =  A  without  loss  of  generality.  Denote  by  P(t)  the  solution 
of  equations  (6.2.1)  and  (6.2.2)  corresponding  to  this  set  of  eigenvalues. 

Let  xi  and  y\  be  n-dimensional  vectors  that  constitute  the  eigvector  corresponding 
to  Ai,  namely, 


Also  for  k  =  0, 1, 2 


$(T,  0) 

V\ 

yi 

=  A 

Xi 

Xl 

(6.4.16) 


.  Let  'y(k)  and  5(k)  be  n-dimensional  vectors  such  that 


7  (k) 

Vi 

.  S{k)  . 

Xi 

(6.4.17) 


Furthermore,  Let  v(k )  = 
we  obtain 


7  (k) 
5(k) 


then  from  equations  (6-4-16)  and  (6.4-17), 


v{k)  =  ${T,0)v(k  +  l)  (6.4.18) 

Define  a  scalar  w(k)  by 

w(k)  =  d*(k  +  1)7 (k  +  1)  -  6*(k)'y(k)  (6.4.19) 

since  ^f(k),  8(k)  — >  0  as  k  — >  00,  we  have 

OO 

Y^w(k)  = -xlyi.  (6.4.20) 

k=0 

On  the  other  hand,  similar  steps  used  to  derive  equation  (6.3.38)  lead  to 

w(h)  >  0,  *  =  0,1,2, .  (6.4.21) 
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We  consider  two  cases: 


(1)  t  7^  iT ,  differentiate  V(t)  along  the  trajectory  of  system  (6.2.12), 

V  (t)  =  x(t)'  P(t)x(t)  +  x  (t)'P(t)x(t)  +  x(t)'  P(t)x(t) 

=  x(t)'{P(t)  +  Ac(t)'P(t)  +  P(t)Ac{t))x(t) 

=  -x(t)'H{t)'H(t)x(t)  (6.4.24) 

<  0. 

Hence  inequality  holds  for  all  H,  £2  7^  iT. 

(2)  t  =  iT,  i  =  0, 1, 2, . .  In  the  rest  of  proof,  we  use  i  to  replace  iT  for 

the  sake  of  brevity.  Since  V(t)  =  —x{t)'H(t)'H{t)x(t)  <  0, 

V(i )  =  V((i  —  1)+)  —  f  x(t)' H  (t)' H  (t)x(t)(1t 

Ji—  1 

V(i+ )  =  x(i)'P(i+)x(i ) 

=  x[i)'P(i)x{i)  +  x(i)'QtmP{i)x(i) 

=  V(t)  +x(i)'QtmP(i)x(i).  (6.4.25) 

Hence, 

H(i+)  -V((i-  1)+)  =  x(i)'Qtmp(i)x(i )  -  /  £(T)'j7(T)^H(r):r(T)dT.(6.4.26) 

Ji— 1 

From  equation  (6.2.12),  we  can  explicitly  solve  x(t )  starting  at  arbitrary  x(i  —  1), 
a;(t)  =  $Ac{t,i  -  l)x(i  -  1),  Vt  G  [(i  -  1),  i].  (6.4.27) 

Replace  x(t)  in  equation  (6.4.26)  by  (6.4.27), 
l/(i+)  -  V((i  -  1)+) 

x(i  1)  [$Ac{i,  i  —  1  )Qtmp(i)$ i  1) 

-/  &Ac(T,i  ~  l)H(T)'H(T)®Ac(r,i  -  l)dr]x(i  -  1).  (6.4.28) 

Ji— 1 
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Since  Ac{t) , H (t)' H (t)  and  Qtmp(i)  are  periodic,  §Ac{i,i  —  1)  =  $^C(T, 0),  and 

Qtmp{T )  =  Qtmp{T)  for  all  i  =  0,1,2, . .  The  equation  (6.4.28)  can  be 

further  simplified  as 

V(i*)  -  V((i  -  1)+)  =  a:(i-l)'[^(T.0)Olmp(T)^(T,0) 

-  /  %'aSt<  0)H(t)'H (r)$Je(r,  (l)dr]j:(i  -  1) 

J  o 

<  0,  by  inequality  (6.2.15).  (6.4.29) 


(Necessity): 

Since  the  closed-loop  system  Ac(t)  is  asymptotically  stable  and  P(t)  is  positive 
definite,  the  inequality  (6.2.15)  follows  immediately  from  Propositions  4.3.1  and 
4.3.2.  □ 


6.4.3  Proof  of  theorem  6.2.1 

We  use  the  standard  completion  of  squares  method  while  accounting  for  the  pos¬ 
sible  jumps.  Differentiating  x'(t)P(t)x(t)  along  the  trajectory  of  system  (6.1.3), 
we  obtain 

^-x'Px  =  x'(t)P(t)x(t)  +  x'(t)P(t)x(t)  +  x'(t)P(t)x(t) 

UjL 

=  x'{t)[A'P{t)  +  P(t)  +  P{t)A\x(t)  +  2  <  u,  B[Px  >  .(6.4.30) 

Replace  term  [A'P(t)  +  P(t)  +  P(t)A]  by  equation  (6.2.1),  and  use  assumption 
ii), 

~^X'PX  —  -x'(t)C[Cix(t)  +  x'(t)P(t)BiB[P(t)x(t)  +  2  <  u,  B[Px  > 

=  -\\z\\2  +  \\u  +  B[Px\\2.  (6.4.31) 
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Integrating  equation  (6.4.31)  from  [iT,  ( i  +  1)T]  for  some  i,  the  right-hand  side 
(RHS)  of  equation  (6.4.31)  becomes 

RHS  =  -|Mlfir,(i+i)T]  +  \\u  +  BiPx\\[iT,(i+i)T]-  (6.4.32) 

where  the  left-hand  side  (LHS)  of  equation  (6.4.30)  becomes 

LHS  =  x'(iT+)P{iT+)x{iT+)-x'{iT)P(iT)x{iT) 

=  [x(»T)  +  B2v(i)]'P{iT+)[x(iT)  +  B2v(i)] 

- x'{iT)P{iT)x(iT ) 

=  v'(i)B'2 P(iT+)B2v(i)  +  2  <  v(i),  B'2P{iT+)x{iT)  > 
+x(iT)[P(iT+ )  -  P(iT)]x(iT) 

=  -72|bWH2  +  ||(7  2I  +  B2P(iT+)B2)R 2 

x[v(i)  +  (72/  +  B'2P(iT+)B2ylB'2P(iT+)x{iT)}\\2.  (6.4.33) 

Since  the  closed-loop  system  is  stable,  hence  x(t)  ->  0  as  t  -»•  oo.  Assume  the 
initial  condition  x(0)  =  0.  Integrate  above  equation  from  0  to  oo,  we  obtain  the 
following  equation, 

OO 

-72|M|22  +  E  ||(7 2I  +  B2P(iT+)B2y /2 

i=0 

x[v(»)  +  {i2I  +  B'2P{iT+)B2)-lB'2P{iT+)x{iT)]\\2 

=  -|kll!2  +  ll«  +  ^Px|||2.  (6.4.34) 

Taking  a  feedback  controller  u(t)  =  -B[P(t)x(t)  and  since  the  summation  term 
is  always  nonnegative,  we  have  the  following  inequality 

NIL<72|M£.  (6.4.35) 

It  is  easy  to  see  from  (6.4.34)  that  the  worst  case  impulsive  disturbances  is 

v(i)  =  -(72/  +  B’2P(k3+)B2)-1B'2P(iT+)x(iT).  (6.4.36) 
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□ 


The  procedure  to  compute  such  a  feedback  controller  is  outlined  as  follows: 

1) .  Given  7  >  0,  check  if  the  matrix  4>(T,  0)  has  no  eigenvalues  on  the  unit 
circle.  If  not,  increase  7  until  the  condition  is  satisfied. 

2) .  Compute  the  periodic  solution  P(t )  according  to  equations  (6.3.9),  (6.3.13) 
and  (6.3.24). 

3) .  Check  if  the  inequality  (6.2.15)  is  satisfied.  If  not,  go  back  to  step  1)  and 
increase  7  until  this  condition  is  satisfied. 

4) .  Compute  the  feedback  controller  u(t)  according  to  equation  (6.2.3). 

6.5  Implementation  By  A  Digital  Controller 

In  this  section  we  will  address  the  optimal  Hqo  control  for  impulsive  disturbances 
by  using  a  sampled-data  controller.  This  consideration  is  partially  motivated  by 
the  need  for  digital  implementation  of  any  controller.  In  chapter  4,  we  have 
deduced  that  the  (sub)  optimal  solution  and  the  controller  is  time-varying  in 
general.  When  using  a  sampled-data  controller  we  are  able  to  convert  the  con¬ 
trol  problem  into  a  discrete-time  control  problem  under  a  proper  sampling 
condition,  hence  the  control  structure  is  simpler  than  the  continuous  one. 

6.5.1  A  sampled-data  controller 

To  facilitate  exposition,  the  problem  we  consider  here  is  a  simple  case  which 
captures  all  the  essential  features  of  the  general  problem.  The  linear  system  is 
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given  by, 


Figure  6.1:  A  sampled-data  controller 


x(t)  =  Ax(t)  +  Biu(t)  +  ^2  B2vS(t  —  kT\) 

k=0 

z(t )  =  Cix(t)  +  Diu(t) 

y(t)  =  C2x(t)  (6.5.1) 

We  shall  study  the  problem  of  controlling  (6.5.1)  using  a  sampled-data  controller 
k,  and  k  has  the  form  of  Figure  6.1,  where  the  S  is  sampling  operator  with 

sampling  period  T2,  k  is  a  discrete-time  LTI  system.  H  is  the  hold  operator 

(zero-order)  [73,  74,  75].  Suppose  that  the  sampling  period  is  T2.  The  structure 
of  this  system  with  sampled-data  controller  is  similar  to  so  called  multirate 
sampling  systems.  [76]  has  shown  that  given  two  sampling  periods  Tx  and  T2, 
if  the  ratio  of  two  periods  is  a  rational  number  T\/T2  =  N-l/N2,  where  Ni  and 
N2  have  no  common  factor.  Then  there  exists  a  smallest  integer  N  and  a  real 
number  T  such  that  T\  =  TNi/N,  T2  =  TN2/N,  and  N  —  N\N2.  If  the  samples 
are  synchronized,  it  follows  that  the  control  and  measurement  signals  will  be 
constant  over  sampling  periods  of  length  T/N.  Sampling  with  that  period  gives 
a  discrete-time  system  that  is  periodic  with  period  T.  The  system  can  then  be 
described  as  time-invariant  discrete-time  system  if  the  values  of  system  variables 
are  considered  only  at  integer  multiples  of  T.  The  ordinary  discrete-time  theory 
and  state-space  form  can  then  be  applied.  In  our  problem,  we  assume  that 
T\/T2  =  n,  where  n  >  1  is  an  integer,  hence  the  above  assumptions  are  satisfied. 
For  simplicity,  we  assume  T\  =  T2  =  T  in  this  chapter.  The  state  and  observation 
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Figure  6.2:  System  diagram 


z(t)  =  0  1  x(t)  =  Cix(t) 

y(t)  =  i  o  x(t)  =  C2x(t) 


(6.5.28) 


Let  us  apply  the  state-space  formula  to  solve  this  control  problem.  We  consider 
the  output-feedback  case.  Without  loss  of  generality,  we  choose  k,c,m 2,  and  bx 
such  that  the  matrices  are  given  as, 

01  0  r  1  r  i 

A  =  ,Bx  =  B  2=  ,Ci=  01  ,c2=  1  0  .(6.5.29) 

-2-3  1  L  J  L  J 

It  is  easy  to  check  that  assumptions  A2 )  and  v43)  are  satisfied.  Let  us  follow  the 
design  procedures. 

Step  1:  Let  sampling  period  T  =  1,  the  M,  and  Pi  are  given  by, 

0.177113  -0.0540768  -0.0885563 

Mi=  ,  Pi  = 

-0.0540768  0.147066  0.0270384 


M2  =  0.127396,  P2  = 


-0.0540768 


0.147066 
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M3  =  0.147066,  P3  =  0.0270384. 


(6.5.30) 


Step  2:  Let  7  =  1,  7 2I  —  M3  =  0.852934  >  0,  and  Q  is  given  by 


Q  = 


0.117005  -0.0410887 

-0.0410887  0.164588 


>0, 


Step  3:  Calculate  A,  Bi,B2,C2  according  to  (6.5.16). 

Step  4:  Solve  two  algebraic  Riccati  equations  (PI)  and  (R2),  the  solutions  are 


(6.5.31) 


It  is  easy  to  show  that  both  M  and  N  are  positive  definite  and  p(NM )  <  1. 
Hence,  the  discrete-time  observer  x(k )  and  suboptimal  controller  u(k)  are  given 
by, 


0.2840 

-0.0086 

0.0405 

-0.0007 

M  — 

-0.0086 

0.1789 

,  N  = 

-0.0007 

1.1872 

x(k  +  1)  = 


0.456713  0.302102 

-0.696961  0.142489 


x(k)  + 


0.145365 

0.0319954 


(y(k)  - 


1  1 


x(k))  + 


u(k), 


u(k)  = 


0.684444  -0.382274 


0.199788 
0.264245 
x(t),  k  =  0,1,2, 


(6.5.32) 


6.6  Conclusion 

In  this  chapter,  we  described  a  complete  state-space  solution  to  the  (sub)optimal 
control  problem  of  a  class  of  linear  systems  subject  to  impulsive  disturbances. 
The  state  feedback  controller  can  be  computed  in  terms  of  the  unique  positive 
definite  periodic  solution  of  a  coupled  Riccati  equations.  Finally,  we  implemented 
this  control  algorithm  by  using  a  sampled-data  controller,  we  showed  that  the 
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structure  of  the  sampled-data  controller  could  be  simpler  that  continuous  version 
under  a  proper  choice  of  the  sampling  period. 
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Chapter  7 


CONCLUSIONS  AND  FUTURE  WORK 


In  this  dissertation,  we  have  derived  a  general  nonlinear  model  for  impact  dy¬ 
namics  on  flexible  structures  based  on  Hertz  law  of  impact.  A  mathematical 
analysis  of  this  model  has  been  carried  out  to  prove  the  existence  of  a  unique 
solution.  A  numerical  method  has  been  developed  based  on  the  contraction 
mapping  principle.  By  taking  advantage  of  the  fact  that  the  impact  period  is 
usually  very  short  we  have  developed  a  series  of  approximation  solutions.  The 
first  order  approximation  yields  a  special  function  which  can  be  used  for  ana¬ 
lytical  and  computational  purposes.  The  second  order  approximation  leads  to 
a  two-parameter  family  of  ordinary  differential  equations  of  which  the  solutions 
contain  universal  features  of  impact  problems.  Simulation  results  of  various  ex¬ 
amples  have  demonstrated  the  usefulness  of  the  developed  numerical  method  and 
approximate  methods.  We  believe  this  nonlinear  model  and  its  approximation 
solutions  are  useful  in  the  design  and  analysis  of  various  engineering  problems 
involving  impact.  This  methodology  has  already  been  used  for  the  analysis  and 
design  of  a  smart  motor  and  a  walking  robot  at  the  University  of  Maryland  at 
college  park.  Besides  the  applications  discussed  in  the  introduction  to  this  dis¬ 
sertation,  the  potential  areas  could  be  found  in  biomechanical  engineering  such 
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as  the  need  to  study  of  human  jumping,  foot  landing  etc. 


We  have  systematically  investigated  the  control  of  impact  forces  by  using  a 
nonlinear  model  for  impact  developed  in  the  first  part  of  this  dissertation,  where 
the  impact  forces  are  treated  as  disturbances  to  the  system.  A  nonlinear  optimal 
control  problem  is  formulated  under  the  assumption  that  the  information  of 
impact  forces  is  known  a  priori.  Optimal  control  strategies  were  derived  via  the 
use  of  dynamic  game  theory.  A  more  interesting  problem  is  to  derive  optimal 
solutions  where  sensors  information  on  impacting  velocities  is  not  available  a 
priori.  This  requirement  leads  to  a  nonlinear  H control  problem  which  needs 
further  investigation.  Again,  by  taking  advantage  of  the  fact  that  impact  period 
is  very  short  in  general  we  have  developed  a  series  of  approximate  solutions  for 
the  nonlinear  optimal  control  problem.  We  have  shown  that  the  higher  order 
terms  are  negligible  in  some  cases,  hence  the  approximation  leads  to  a  linear 
problem.  A  linear  H0 0  control  problem  is  formulated  and  a  state-space  solution 
has  been  obtained.  The  solution  is  naturally  associated  with  the  existence  of  a 
stabilizing  periodic  solution  of  a  coupled  Riccati  equations.  Hamiltonian  theory 
is  employed  to  analyze  the  coupled  Riccati  equations.  Finally,  we  investigate 
the  digital  implementation  of  this  control  algorithm  by  using  a  sampled-data 
controller.  We  show  that  under  a  certain  sampling  condition,  the  controller 
structure  could  become  simpler  than  the  continuous  one.  These  control  methods 
are  proposed  and  analyzed  on  the  theoretical  base.  It  would  be  more  interesting 
and  challenging  to  implement  and  test  these  algorithms  on  some  real  world 
problems.  One  application  is  force  and  constrained  motion  control  of  a  flexible 
robot.  An  experimental  apparatus  in  the  ISL  laboratory  at  the  University  of 
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Maryland  could  be  a  good  test  bed.  Another  application  is  to  investigate  the 
active  control  of  the  suspension  system  of  a  vehicle  [16,  15]. 
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Appendix  A 


PROOFS  OF  CHAPTER  3 


A  .1  A  simply  supported  beam 


For  the  free  vibration,  the  deflection  of  a  beam  is  governed  by  the  PDE, 
d2w  d 4w 


dt 2  +  dx4 


=  0 


0  <  x  <  l, 


(A.1) 


and  boundary  conditions  for  the  simply  supported  case  are 


w(0,t)  =  w"(0  ,t)  =  0,  (A. 2) 

w(l,t)  —  w"(l,t )  =  0. 

The  corresponding  eigenvalue  problem  can  be  written  as 

W""  =  0W(x),  {^L  =  W"")  0  <  x  <  l,  (A. 3) 

where  (3  is  an  eigenvalue  and  the  W(x)  €  L2[0,l]  is  the  corresponding  eigen¬ 
function.  The  boundary  conditions  are 

tr(0)  =  w"(0)  =  0,  (A. 4) 

W(l)  =  W"(l)  =  0. 
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The  general  solution  can  be  expressed  as 


W  ( x )  =  c\sin(3x  +  C2Cos/3x  +  cisinhfdx  +  c^coshfdx,  (A. 5) 

where  the  C\  to  C4  are  constants  to  be  determined  by  boundary  conditions  (A. 4). 
In  order  to  satisfy  boundary  conditions  (A. 4)  we  get, 


sin(3l  =  0.  (A. 6) 

The  solution  of  equation  (A.  6)  is  (3kl  =  kn,  k  —  1, 2, . 00 

The  normal  eigenfunctions  are 

Wk(x)  —  \[2sin(3kx  k  =  1,2, . 00  (A. 7) 


Lemma  A  .1  V:c,  (  G  [0,/]  and  Vt  >  0  ,  the  Green’s  function  is  uniformly 
bounded. 

Proof.  It  is  easily  checked  that  (A. 3),  (A. 4)  is  self  -adjoint.  Hence,  its 
Green’s  function  can  be  expressed  as 

G(x,Ct)  =  ^Wt(x)Wt(0S^^H(t).  (A.8) 

fc=l  Wk 

Since, 

Wk(x)  —  V2sin(3kx  Vx  €  [0,  l] 

|W/’fc(^)W/fc(C)|  =  2\sin/3kxsinpk(\  <  2 

k=  1  Wk 

k= 1  Wk 

^  2  \sinwkt\ 

hx  wk 

00  1 

<  2  E  —  («£  =  ftp/EI) 

k=lWk 
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converges,  3 M  >  0,  such  that  G(x,(]t)  <  M  , 


Va:,C  €  [0,  /]. 


A  .2  A  cantilevered  beam 


The  PDE  is  same  as  in  equation  (A.l),  and  the  boundary  conditions  for  the 
cantilevered  beam  are 

w(0,t)  =  w  (0,t)  =  0,  (A. 9) 

w(l,t )"  =  )  =  0. 

The  approach  to  solve  the  eigenvalue  problem  is  the  same  as  in  A.l.  Again,  we 
write  the  general  solution  and  plug  in  the  boundary  conditions  (A. 9)  into  this 
equation.  After  simplification,  we  get  the  beam  characteristic  equation, 


cos/3lcosh/3l  —  —  1. 


(A.10) 


The  orthonormal  eigenfunctions  are  determined  by  the  following  equations, 


Wk(x) 
where  Wk(x) 
and  A\ 


i^x) 

cosh/3kx  —  cos(3kx  sinh/3kx  —  sin(3kx 
cosh/3kl  +  cos(3kl  sinhf3kl  +  sinf3kl 


(A.ll) 


=  /  Wjj(x)dx 

JO 


cos2f3kl 

sin40kl 


*  =  1,2, 


oo. 


Lemma  A  .2  The  orthonormal  eigenfunctions  { Wk{x )}£Li  are  uniformly 
bounded. 


Proof:  There  are  infinitely  many  solutions  to  the  characteristic  equation  (A.  10), 
0  <  flil  <  fol  < . <  oo,  where  (5\l  —  4.73. 

Note  that  cosh/3kl  >  0,  coshfikx  >  0,  sinhpkl  >  0,  sinh(3kx  >  0.  Vx  6  [0,  Z],VA;. 

W  ^  _  cosh(3kx  —  cos(3kx  sinh(3kx  —  sin(dkx 

coshf3kl  +  cos(3kl  sinhfikl  +  sin(3kl 
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_  _ Ck(x)  -  Dk (x) _ 

( cosh(3kl  +  cos/3kl)(sinhfikl  +  sin(3kl ) 
where  Ck(x)  =  (cosh/3kx  —  cosfikx)(sinhfikl  +  sinj3kl); 

Dk(x)  =  ( sinh[3kx  —  sinfikx)(cosh(5kl  +  cos(5kl). 

Now,  we  carry  out  the  simplifications: 

Ck(x)  —  Dk(x)  —  coshfikxsinh/3kl  —  sinh/3kxcosh/3kl  —  cosfikxsinh/3kl 
+cosh/3kxsin/3kl  +  sin/3kxcosh(3kl  —  sinh/3kxcos/3kl 
+sin(5kxcosfikl  —  cos(3kxsin@kl 
=  sinhfik(l  —  x)  +  sin/3k(x  —  l)  —  cosfikxsinhfikl 
+coshfikxsinfikl  +  sinj3kxcoshfikl  —  sinhfikxcosfikl 


|  Ck(x)  —  Dk(x)  |  <  sinhfik(l  —  x)  +  \sin/3k(x  —  Z)|  +  \cos/3kx\sinh/3kl 

-\-cosh/3kx\sinf3kl\  +  \sin(3kx\cosh(3kl  +  sinhfikx\cosfikl\ 

<  sinhfik(l  -  x)  +  1  +  sinhfikl 
+cosh/3kx  +  coshfikl  +  sinhfikx 

<  l/2ew~x)  +  1  +  ePkl  +  ePkX  =  hk(x)  >  0; 


Since  (3kl  >  4  \/k,  it  follows  that  sinh(5kl  >  2,  and  coshfikl  >  2 ,Mk.  Hence, 


\Wk\  < 

< 

< 

< 


\Ck(x)  -  Dk(x) 


\coshfikl  +  cos(3kl\\sinhfikl  +  sin/3kl\ 

_ hk(x) _ 

(cosh/3kl  —  1  )(sinh/3kl  —  1) 

2  hk(x) 


coshj3kl(sinh/3kl  —  1) 
4  hk(x) 

cosh/3klsinh(3kl 


since  coshfikl  >  2; 


since  sinh{3kl  >  2; 


since  cos/3klcoshf3kl 


-1 


cos2  /3kl 


i 

cosh2(3kl 


and  cosh(3kl  >  0; 
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For  the  cantilevered  beam,  the  differential  operator  is  also  self-adjoint.  Thus, 
the  proof  that  the  Green’s  function  is  uniformly  bounded  is  similar  to  lemma 
A.l. 
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